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A  study  of  particle  collision  and  concentration  non-uniformity  in  turbulent  flows 
is  carried  out.  A  correction  to  the  classical  result  of  Smoluchowski  for  the  collision  rate 
of  monodisperse  particles  in  a  laminar  shear  flow  has  been  made.  A  theoretical 
framework  has  been  developed  to  evaluate  the  collision  rate  and  collision  velocity  of  very 
small  particles  in  general  turbulent  flows.  The  present  approach  differs  significantly 
from  the  classical  approach  of  Saffman  &  Turner  in  that  the  ensemble  average  is  taken 
after  the  collision  rate  for  a  given  flow  realization  is  calculated.  This  avoids  the 
assumption  of  isotropy,  as  needed  in  the  classical  approach,  and  allows  for  the  evaluation 
of  the  collision  rate  in  general  turbulent  flows. 

Direct  numerical  simulations  and  analytical  studies  on  particle  collisions  in  an 
isotropic,  Gaussian  turbulence  have  been  carried  out.  Several  different  collision  models 
(or  post-collision  treatments)  are  implemented  in  the  simulations  and  the  difference 
among  various  collision  models  is  examined.  The  effects  of  particle  inertia,  size,  gravity, 
and  mean-concentration  on  the  collision  rate,  collision  velocity,  and  collision  angle  are 
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investigated.  Asymptotic  analyses  have  been  developed  to  predict  the  effects  of  small 
inertia  on  the  increases  in  collision  rate  and  collision  velocity.  Asymptotic  prediction  for 
the  increased  collision  velocity  agrees  well  with  the  results  of  both  numerical  simulation 
and  PDF-based  method.  The  increased  collision  velocity  is  found  to  contain  terms  of 

both  0((Pxk)  )andO((Pxk)  ),  where  p  is  the  particle  response  time  and  xk  is  the 
kolmogorov  time  scale.  The  increases  in  particle  size,  gravitational  acceleration,  and 
mean-concentration  in  general  reduce  the  normalized  particle  collision  rate  in  isotropic 
turbulence  for  monosized  particles  as  opposed  to  very  small  particles  in  the  very  dilute 
limit. 

A  new  method  has  been  developed  to  quantify  the  concentration  non-uniformity 
of  particles  in  turbulent  flows.  The  effects  of  particle  size  and  inertia  on  the  concentration 
non-uniformity  are  examined  and  elucidated.  The  concentration  non-uniformity  can  be 
greatly  overestimated  if  the  effect  of  finite  particle  size  is  neglected. 
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CHAPTER  1 
INTRODUCTION 


Particles  and  turbulence  are  around  us  everywhere  and  greatly  impact  our  daily 
life.  Dusts,  sand,  rain  drops,  fine  powders,  glass  beads,  viruses  and  bacteria  are  all 
'particles.'  Most  of  the  fluid  flows  occurring  in  nature  such  as  air  current  and  water  flows 
in  rivers  and  seas  are  turbulent  flows.  The  particles  in  the  atmosphere  remain  airborne  for 
long  periods,  lowering  visibility.  The  brown  haze  that  is  often  seen  over  large  cities  in 
autumn,  winter  and  spring  is  due  mainly  to  the  interactions  between  particles  and  the  air. 

The  improvement  of  particulate-fluid  systems  requires  a  better  understanding  of 
the  fundamental  physics  of  the  particle-particle  and  particle-fluid  interactions.  The 
behavior  of  particulate-fluid  mixtures  may  vary  widely  depending  upon  the  conditions 
under  which  the  particle-fluid  system  functions.  Particle-particle  interactions  involve  the 
collisions  among  particles  which  depend  on  inter-particle  forces,  particle  size  distribution, 
and  the  concentration  distribution.  Particle-fluid  interactions  include  the  hydrodynamic 
forces  acting  on  particles  by  fluid  and  the  feedback  to  fluid  from  particles,  etc.  Particle 
collision,  coagulation,  and  aggregation  are  the  examples  of  most  important  scenarios  of 
those  interactions  that  are  determined  by  particle  properties  and  flow  structures.  The  fact 
that  collisions  do  occur  among  the  particles  of  size  ranging  from  a  few  Angstroms  to  a 
few  millimeters  explains  the  vast  amount  of  scientific  and  industrial  applications  of  the 
collision  process.  For  example,  molecular  collisions  help  explain  the  physical  properties 
such  as  viscosity,  diffusivity,  and  conductivity  in  gases  and  liquids.  Numerous  industrial 
applications  can  be  cited:  design  of  fine  spray  combustion  nozzles,  control  of  industrial 
emissions  (Presser  &  Gupta  1993,  Tu  et  al.  1996)  and  the  understanding  of  the  collision 
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process  among  suspended  particles  in  pneumatic  transport  systems  or  in  gas  cleaning 
chambers  (Abrahamson  1975,  William  &  Crane  1983),  and  in  sewage  disposal  devices 
(Hunt  1982),  where  particle  size  distribution  and  concentration  non-uniformity  are 
important  characteristics  that  depend  on  particle  collision  and  coagulation.  Most  of 
previous  studies  on  particulate  removal  have  focused  on  the  filters,  with  little  attention  to 
the  particle-particle  interactions  and  particle-fluid  interactions.  Despite  many  years  of 
research  and  development  in  particulate  removal  technologies,  the  collection  of  small 
particles  with  size  less  than  1  um  is  still  a  rather  difficult  task.  It  is  extremely  difficult  to 
remove  fine  organic  active  particles  (<0.03  um)  from  a  waste  stream  because  the  particles 
are  much  smaller  than  the  pores  in  existing  filters.  Typical  collection  efficiencies  of  the 
filtration  devices  decrease  rapidly  with  the  decrease  of  particle  size  in  the  range  0.2-30um 
(Davies  1953).  A  good  alternative  is  to  aggregate  fine  particles  which  can  be  removed 
using  conventional  technologies.  For  the  coagulation  and  aggregation  to  occur,  particles 
have  to  collide  with  each  other  and  adhere.  It  is  very  desirable  that  one  can  predict  how 
the  small  particles  are  driven  to  collide  with  each  other  so  they  aggregate  or  grow  in  size 
and  become  large  enough  to  be  captured  by  the  filters.  Obviously,  the  efficiency  of 
particle  removal  mainly  depends  on  the  control  of  particle  collisions  and  concentration 
non-uniformity  in  filtration  systems  (Fathikalajahi  et  al.  1996).  It  is  evident  that  particle 
collision  is  one  of  the  most  important  interactions  for  understanding  the  fundamentals  of 
particle-fluid  systems. 

Particle  collision  kernel  and  collision  rate  are  critical  quantities  in  the 
development  of  collision  theories.  The  rate  of  flow  induced  particle  collisions  affects 
directly  the  particle  size  distribution  in  particle  production.  Examples  of  turbulence 
induced  particle  collision  can  be  seen  in  the  formation  of  rain  drops  in  clouds,  pulverized 
coal  combustion,  agglomeration  of  fine  powders  in  gas  flows,  air  filtration  equipment, 
sewage  disposal  devices,  fast  fluidized  beds,  dust  and  spray  burners,  and  so  on.  Particle 
collisions  can  take  place  through  a  variety  of  collision  mechanisms.  Particle  collisions 
may  be  induced  by  centrifugal  and  gravitational  forces,  thermal  forces  in  Brownian 
motion,  and  van  der  Waal  forces.  In  an  actual  system,  the  scenario  after  every  contact 
between  particles  may  be  different  from  one  to  another,  depending  on  the  underlying 
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physics  of  the  system.  The  collisions  among  raindrops  enhance  the  growth  of  droplets  in 
the  air,  since  the  contact  between  droplets  may  most  likely  result  in  the  coalescence  of 
two  droplets.  The  growth  of  bubbles  in  a  liquid  flow  is  another  example  of  particle 
coagulation  upon  collisions.  The  contacts  between  two  solid  particles  may  behave  in 
various  different  ways.  When  two  solid  particles  are  brought  together  the  colliding  pair 
may  stick  or  separate  after  collision  upon  contact.  Many  theoretical  efforts  were  made  to 
predict  the  collision  kernel  and  rate.  To  validate  a  theoretical  analysis,  direct  numerical 
simulations  were  carried  out  to  obtain  statistical  average  of  the  collision  rate  in  a  particle- 
fluid  system.  Different  interpretations  of  the  immediate  consequence  after  the  contact 
between  a  pair  of  particles  have  resulted  in  different  collision  models  or  so-called  "post- 
collision  treatment"  schemes.  To  quantify  the  post-collision  treatment,  a  statistical 
parameter  called  collision  efficiency  has  been  introduced.  The  collision  efficiency  is  one 
if  particles  coagulate  to  form  a  larger  one.  Usually,  it  is  less  than  one  since  some  of  the 
colliding  particles  may  bounce  back.  Although  many  theoretical  and  experimental  studies 
have  been  devoted  to  estimates  of  the  collision  efficiency  for  water  droplets  settling 
through  still  air  (Klett  &  Davis  1973,  Lin  &  Lee  1975,  Beard  &  Ochs  1983)  and  droplets 
settling  down  through  laminar  shear  flow  (Jonas  &  Goldsmith  1972),  even  an 
approximate  estimate  of  the  collision  efficiency  in  a  complicated  situation  can  be  very 
difficult  to  obtain.  At  this  stage,  it  is  still  not  possible  to  obtain  an  accurate  solution  for 
the  evaluation  of  collision  efficiency  in  complex  systems.  To  better  focus  on  the  issues  of 
our  interests,  the  collision  rate,  the  collision  efficiency  is  assumed  to  be  either  one  or 
zero  in  this  study.  Although  a  unity  collision  efficiency  will  overestimate  the  size  growth, 
however,  theoretical  and  numerical  exploration  on  the  governing  physics  may  still  hold 
its  significance  for  a  better  understanding  of  the  interactions  in  particle-fluid  systems. 

Population  balance  equations  are  often  used  for  predicting  the  evolution  of 
particle  size  distribution  in  particle  production  systems  or  in  studying  particle 
coagulation.  They  are  typically  given  in  the  form  of 
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dnj 

"dT  =  -aUn\n\  -  a12nln2  "  a13nln3  '  - 


dn2  i 

dt  =  2  ^  1n1n1  -  a12n,n2  -  a22n2n2-  a23n2n3  - ... 


(1.1) 


dn3 

"dT  =  a12nln2  "  a13nln3  -  a23n2n3  "  <*33n3n3  "  - 


dn4  1 

dt    ~  a13nln3  +  2  a22n2n2  *  a14nln4  "  a24n2n4  "  - 


In  the  above,  n;  is  the  number  concentration  of  the  z'th  particle  size  group.  For  primary 
particles,  i=l,  with  a  volume  vj,  the  particle  in  the  rth  group  has  a  volume  of  ivx.  The 


is  called  collision  rate  which  represents  the  number  of  collisions  among  rth  particle  with 
yth  particle  per  unit  volume  and  per  unit  time.  The  coefficient  is  called  the  collision 
kernel  or  collision  function,  it  depends  on  particle  size,  particle  inertia  and  flow  structure, 
and  must  be  evaluated  separately.  The  collision  rate  or  collision  kernel  is  a  measure  of 
the  ability  of  particles  to  collide,  coagulate,  and  aggregate.  It  is  essential  to  determine  the 
instantaneous  particle  size  distribution  and  momentum  transfer  among  particles  in  a 
particulate  system.  Over  past  80  years,  many  efforts  have  been  devoted  to  developing  the 
particle  collision  theories  for  predicting  the  particle  collision  rate  in  fluid  flows.  A  brief 
review  is  given  below. 


A  geometrical  collision  occurs  when  two  particles  reach  a  separation  distance  less 
than  R=rj+rj,  in  which  r\  and  rj  are  the  respective  radii  of  two  particle  groups  and  R  is  the 
radius  of  the  colliding  surface.  Theoretical  analysis  of  geometrical  collisions  in  a  uniform 
shear  flow  was  first  carried  out  by  Smoluchowski  (1917).  At  a  given  instant  the  particle 


quantity  dyn^,  denoted  as  nCy  in  this  thesis, 


nCij  -  aynjiij 


(1.2) 


1 . 1  Smoluchowski's  Prediction  in  Laminar  Shear  Flow 
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number  concentrations  per  unit  volume  are  nj  and  nj  for  each  size  group.  Taking  any 
particle  of  the  fth-group  as  a  target  particle,  the  volume  flow  rate  of  the y'th-group  particles 
reaching  the  colliding  surface  of  radius  R  of  the  given  target  particle  is  nj}Wr£0wrdS  in 

which  w  is  the  relative  velocity  between  z'th-group  and  y'th-group  particles  and  wr  is  the 
radial  component  of  w.  Thus  the  number  of  collisions  occurring  among  z'th-group  and 
yth-group  particles  per  unit  time  per  unit  volume  is  given  by 

n  cij  =  njnj  JWf  s0  wrdS  =  n^j  j  r(ri+rj)3  =  ay  njnj  (1.3) 

where  T  is  the  shear  rate  and 

air  yHri+rj)3  (1.4) 

is  often  called  the  collision  kernel.  For  like  particles,  the  number  of  collisions  per  unit 
time  and  per  unit  volume  was  often  simply  taken  as 

no^nin^r^)3.  (1.5) 
The  factor  ^  was  introduced  to  avoid  double  counting  of  the  collisions  among  the  like 
particles.    When  the  above  results  are  used  in  the  population  balance  equations,  for 

example,  the  rate  of  the  "destruction"  of  z'th-group  particles  due  to  collisions  is  -2nCii  + 

4        3  4  3 

I  n  c  =  -  njnj  -  T(2rj)  -  I  njnj  -  r(rj+rj)J  since  each  collision  among  the  like  particles 

j*i  3  j*i  3 

results  in  the  loss  of  two  particles  of  that  size  group.  Equation  (1.5)  has  been  extensively 
used  in  the  literature  in  calculating  the  particle  collision  rate  and  particle  size  evolution. 

In  addition  to  Smoluchowski's  work,  Hocking  &  Jonas  (1970)  also  calculated  the 
collision  rate  for  the  sedimenting  particles  in  a  uniform  shear  flow  using  Stokesian 
theory.  They  studied  the  collision  efficiency  by  which  the  particle  coagulation  was 
assumed  if  the  surface  distance  between  two  neighboring  particles  was  less  than  a 
prescribed  fraction  of  particle  diameter. 

1.2  Various  Theories  for  Collision  Rate  in  Turbulence 


When  particles  collide  in  complex  flows,  such  as  turbulent  flows,  or  when  more 
than  one  collision  mechanism  are  involved  in  the  particulate  system,  it  is  quite  difficult  to 
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derive  an  exact  expression  for  collision  rate.  Numerous  efforts  have  been  made  to 
estimate  the  collision  rate  by  means  of  evaluating  the  collision  velocity  in  turbulence 
during  past  decades.  Camp  and  Stein  (1943)  employed  the  concept  of  mean  velocity 
gradient  related  to  the  turbulent  energy  dissipation  rate  e  and  kinematics  viscosity  v  and 
obtained  the  collision  kernel  for  turbulence  shear-induced  collision  rate  based  on 
Smoluchowski  work  (1917)  in  laminar  shear  flow. 

a„=^R3A1/2.  (1.6) 
IJ     3  v 

Saffrnan  and  Turner  (1956)  gave  two  predictions  of  particle  collision  kernel  in  an 
isotropic,  Gaussian  turbulence.  In  the  first  prediction,  it  was  assumed  that  particles  are 
inertialess  so  that  they  follow  the  fluid  completely  and  particle  size  is  small  compared 
with  the  Kolmogorov  length  scale  of  the  turbulence.  The  collision  kernel  was  given 

a  =1294  R3(-)1/2.  (1.7) 
u  v 

Their  second  prediction  was  for  the  particles  with  small  inertia  and  gravity. 


Po  Dt 

+  {(l-^)2(T,-T2)V^R2-],/2.  (1-8) 

3       p0  9  v 

However,  collision  kernel  for  like  particles  or  the  particles  with  same  inertia  is  reduced  to 

a..  =1.67  R3(-),/2.  (1.9) 

Thus  the  effect  of  particle  inertia  and  gravity  disappeared  in  the  monodisperse  case, 
which  is  not  physically  sound.  In  addition,  the  result  of  their  second  analysis  was 
inconsistent  with  the  first  one,  as  easily  seen  from  (1.7)  and  (1.9). 

Abrahamson  (1975)  extended  the  concept  of  kinetic  theory  of  gases  for  the 
collisions  among  large  particles  in  a  high  Reynolds  number  turbulence.  It  was  assumed 
that  the  velocities  of  colliding  particles  are  uncorrelated,  normally  distributed  in  the 
absence  of  gravitational  settling.  The  variance  of  particle  velocity  v  was  estimated  by 
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<v2>=   (1.10) 

1+  15  xp  s/u2 

where  u  is  the  mean-squared  fluctuating  velocity  of  the  fluid  and  xp  is  the  particle 
response  time.  The  collision  kernel,  without  gravitational  settling,  was  given  by 

ay  =  5.0  R2  (  <v2>  +  <v?j>  )1/2.  (1.11) 

Balachandar  (1988)  performed  an  analysis  on  the  collision  kernel  by  assuming 
that  the  probability  distribution  function  (PDF)  of  the  velocity  gradient  in  turbulence 
follows  a  log-normal  PDF  distribution.  He  found  a  correction  to  Saffman  &  Turner's 
prediction  for  collision  kernel  in  isotropic  turbulence.  Williams  &  Crane  (1983)  and 
Kruis  &  Kusters  (1996)  performed  stochastic  analysis  to  obtain  the  collision  kernel  with 
inclusion  of  the  combined  effects  of  local  shear  and  particle  inertia.  Zhou  et  al.  (1997) 
carried  out  a  similar  analysis  and  modified  the  Kruis  &  Kusters'  result  in  turbulence. 
Those  analyses  on  the  collision  kernels  were  obtained  on  the  basis  of  two-point  velocity 
correlations  without  tracking  the  trajectories  of  individual  particles  in  the  turbulent  flow 
(hereinafter  referred  as  Eulerian  method).  However,  we  have  found  that  the  collection  of 
the  relative  velocities  (or  collision  velocities)  upon  collisions  based  on  tracking  the 
trajectories  of  individual  particles  in  turbulent  flow  is  a  biased  average  in  favor  of  those 
regions  where  the  collision  velocity  is  high.  Hereinafter,  this  method  based  on  particle 
tracking  is  referred  as  Lagrangian  method.  The  detailed  discussions  on  the  difference 
between  the  results  of  Eulerian  and  Lagrangian  methods  on  collision  velocity  will  be 
presented  in  chapter  IV. 


1 .3  Recent  Development  in  Computer  Simulation  on  Particle  Collision  Rate 

Recent  developments  of  direct  numerical  simulations  (DNS)  in  turbulence  have 
greatly  enhanced  the  understanding  of  the  turbulence-induced  particle  collisions. 
Balachandar  (1988)  also  carried  out  DNS  of  particle  collision  in  an  isotropic  turbulence. 
It  was  found  that  the  probability  distribution  function  of  the  velocity  gradient  follows  a 
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log-normal  distribution  instead  of  Gaussian  distribution.  The  collision  kernel  was 
obtained  in  a  frozen  turbulence  with  a  maximum  of  500  particles  in  the  flow  field. 
Artificially  large  particles  were  used  in  order  to  increase  the  collision  events  within  a 
reasonable  computational  time.  It  was  concluded  that  the  particle  inertia  does  not  affect 
the  average  collision  kernel  for  large,  monosize  particles.  Large  discrepancy  in  the 
collision  rate  for  monosize  inertialess  particles  between  his  numerical  results  and 
Saffman  &  Turner's  prediction  was  reported.  However  it  is  not  clear  whether  this 
discrepancy  is  due  to  the  difference  in  the  small  scale  turbulence  structure  or  not. 

Chen  et  al.  (1995)  carried  out  DNS  results  for  the  aerosols  in  a  turbulent  channel 
flow.  They  mainly  focused  on  the  possible  effects  of  collisions  on  the  deposition  rate  of 
the  droplets  on  the  channel  walls.  They  found  that  the  particle  concentration  in  the  wall 
region  can  be  significant  while  the  bulk  concentration  is  small.  The  collision  kernel  and 
rate  were  found  to  peak  in  the  wall  region,  and  particle  depositions  can  be  enhanced  by 
the  formation  of  larger  particles. 

Wang  et  al.  (1997)  also  carried  out  DNS  to  investigate  the  collision  kernel  in 
homogeneous  turbulence.  They  allowed  particles  in  the  system  to  overlap  in  space  which 
is  unphysical.  The  effects  of  finite  inertia  on  collision  kernel  were  also  reported  in  their 
DNS  results  (Zhou,  Wexler  &  Wang  1997).  Zhou  et  al.  (1997)  developed  an  analysis 
showing  that  the  increases  in  particle  collision  velocity,  the  non-uniformity  of  particle 
concentration  and  the  compressibility  of  particle  velocity  in  the  continuum  sense  are  three 
possible  mechanisms  to  affect  the  overall  collision  kernel  for  particles  with  finite  inertia. 
However,  they  did  not  quantify  these  terms  separately  in  their  work.  Thus  how  the 
collision  velocity  and  the  concentration  non-uniformity  play  in  enhancing  the  overall 
collision  kernel  remains  unknown  at  this  stage. 


1 .4  Recent  Development  in  Particle  Concentration  Non-uniformity 


Recent  studies  have  shown  that  the  particle  concentration  in  turbulence  may  be 
highly  nonuniform.  Heavy  particles  with  finite  inertia  tend  to  accumulate  in  the  regions 
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of  high  strain  rate  and  to  move  away  from  the  vortex  cores.  The  occurrence  of  particle 
accumulations  is  referred  as  concentration  non-uniformity  or  preferential  concentration. 
The  coherent  vortical  structures  are  the  mechanisms  that  cause  the  nonuniform 
distribution  of  particles.  The  coherent  vortices  and  concentration  non-uniformity  have 
been  observed  by  many  numerical  simulations  and  experimental  visualizations.  The 
numerical  simulations  for  demonstration  of  the  concentration  non-uniformity  were  carried 
out  in  plane  mixing  layers  (Crowe  et  al.  1988,  Chien  &  Chung  1987,  Wen  et  al.  1992, 
Tang  et  al.  1992,  Wang  1992),  in  wake  flows  (Chien  &  Chung  1988,  Tang  et  al.  1992), 
in  jet  flows  (Chung  &  Trout  1988,  Hansell  1992),  in  wall-bounded  flows  (Kallio  &  Reeks 
1989,  McLaughlin  1989,  Yonemura  et  al.  1993,  Rouson  &  Eaton  1994,  Tanaka  et  al. 
1996),  in  homogeneous  flows  (Maxey  &  Corrsin  1986,  Hunt  et  al.  1987,  Fung  &  Perkins 
1989,  Squires  &  Eaton  1991,  Wang  &  Maxey  1993,  Fessler  et  al.  1993).  The  non-uniform 
distribution  of  particles  has  also  been  found  by  the  experiments  in  plane  mixing  layers 
(Kobayashi  et  al.  1988,  Kamalu  et  al.  1988,  Lazaro  &  Lasheras  1992,  Ishima  et  al.  1993), 
in  wake  flows  (Tang  et  al.  1992,  Yang  et  al.  1993,  Bachalo  1993),  in  jet  flows  (Longmire 
&  Eaton  1992),  in  wall-bounded  flows  (Rashidi  et  al.  1990,  Young  &  Hanratty  1991). 
Wang  &  Maxey  (1993)  have,  through  direct  numerical  simulations  (DNS)  in  isotropic 
turbulence,  shown  that  the  preferential  concentration  achieves  maximum  at  Pxk  =1  in 
which  p  is  the  reciprocal  particle  response  time  and  xk  is  the  Kolmogorov  time  scale  of 
the  turbulence.  One  of  the  consequences  of  the  preferential  concentration  is  the 
enhancement  of  the  settling  velocity  of  heavy  particles  in  the  Stokes  drag  range  (Maxey 
1987,  Mei  1994).  One  of  the  related  industrial  applications  is  that  the  particulate 
removal  efficiency  of  the  venturi  type  scrubber  was  affected  by  the  nonuniform 
concentration  distribution  of  droplets  (Fathikalajahi  et  al.  1996). 

Particle  trajectories,  therefore  individual  particle  positions,  can  be  obtained 
accurately  only  in  Lagrangian  frame  by  solving  the  equations  of  particle  motions. 
Particles  have  finite  response  times  in  general  but  they  do  not  have  an  intrinsic  diffusivity 
or  viscosity  associated  with  the  motion  of  individual  particles.  Turbulence  particle 
diffusivity  is  the  consequence  of  ensemble  averaging  over  many  realizations  and  it  is  used 
for  predicting  the  variation  of  average  concentration  field  on  a  large  scale.  It  cannot  be 
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used  to  describe  the  small-scale  spatial  variation  of  the  particle  concentration.  From  the 
computational  point  of  view,  a  DNS  using  1283  grid  resolution  for  an  isotropic 
turbulence  becomes  routine  and  DNS  using  2563  or  even  5 12 3  have  been  carried  out. 
For  a  reasonable  contour  representation  of  particle  concentration,  it  is  desirable  to  have  at 
least,  say,  10  particles  in  each  cell  so  that  the  statistical  noise  does  not  become 
overwhelming.  This  requires  1.7xl08  particles  with  a  2563  grid  resolution.  Due  to  the 
preferential  concentration,  there  will  be  regions  and  cells  that  do  not  contain  any  particle 
or  only  one  or  two  particles  at  a  given  instant.  Thus,  an  accurate,  three-dimensional 
contour  representation  of  instantaneous  particle  concentration  with  such  high  resolution 
will  still  be  heavily  contaminated  by  the  inherent  statistical  noise.  No  such  three- 
dimensional  concentration  contour  has  been  reported  to  date.  Due  to  the  computer  power 
limitation,  previous  numerical  simulations  involved  anywhere  from  104  to  106  particles, 
which  is  far  below  the  level  of  reaching  a  reasonable  statistics. 

It  has  been  difficult  to  quantify  the  extent/degree  of  particle  concentration  non- 
uniformity  induced  by  the  small-scale  turbulence  structure.  Squires  &  Eaton  (1991)  and 
Wang  &  Maxey  (1993)  showed  in  their  DNS  the  similarity  between  the  instantaneous 
particle  concentration  contours  with  the  contours  of  the  vorticity  magnitude  (or 
enstrophy)  in  the  same  plane  at  the  same  instant.  To  quantify  the  deviation  of  the  particle 
concentration  field  from  a  statistically  uniform  distribution,  Wang  &  Maxey  (1993)  used 
a  statistical  quantity  Z(P(k)-  Pb(k))2  averaged  over  space,  where  P(k)  is  the 
probability  of  finding  k  particles  in  a  given  cell  and  Pb(k)  follows  the  Bernoulli 
distribution  for  a  random  distribution  of  the  particles.  While  this  quantity  shows  that  the 
preferential  concentration  maximizes  near  Pxk=l  (Wang  &  Maxey  1993),  it  varies  with 
total  number  of  particles  in  each  of  their  simulations.  The  particle-particle  interactions 
are  neglected  in  their  simulation  and  no  systematic  investigation  on  its  dependence  on  the 
cell  size  was  carried  out.  Hence,  its  usefulness  is  limited.  Although  numerous 
experimental  studies  have  clearly  shown  the  non-uniform  particle  distribution  in  turbulent 
flows,  it  is  quite  difficult  to  measure  the  3D  non-uniformity  so  as  to  isolate  the  effects  of 
flow  shear,  particle  inertia  and  gravity. 
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Obviously,  the  local  accumulation  will  affect  the  encounters  among  the  particles 
and  thus  the  overall  collision  rate  is  expected  to  be  modified  accordingly.  However,  the 
effect  of  the  preferential  concentration  on  the  collision  kernel  has  never  been  estimated  or 
separated  from  other  factors  such  as  the  increased  particle  collision  velocity  due  to  small 
particle  inertia. 


1.5  Scope  of  Present  Study 

Particle  collisions  are  extremely  difficult  to  observe  experimentally,  even  with 
today's  sophisticated  laser  instrumentations.  Direct  numerical  simulation,  however,  has 
its  advantages  of  providing  detailed  information  on  particle  trajectories,  flow  structures 
and  is  ideally  suited  for  parametrical  investigations  of  particle  and  fluid  flow  properties. 
Numerical  implementation  of  particle  interactions  in  a  simulation,  while  posing  a 
considerable  numerical  challenge,  introduces  no  fundamental  difficulties.  In  this  work, 
extensive  data  obtained  by  numerical  simulations  on  the  collision  velocity,  collision 
angle,  collision  kernel  and  collision  rate  are  vital  to  understanding  of  the  fundamental 
physics  in  the  particle-fluid  interactions.  One  of  the  major  objectives  of  this  thesis  is  to 
develop  analytical  and  numerical  methods  for  predicting  the  particle  collision  rate  and 
concentration  non-uniformity  with  respect  to  the  effects  of  turbulent  shear,  particle  inertia 
and  gravity  in  particle-fluid  systems. 

As  mentioned  in  §1.2,  previous  theories  on  the  collision  rate  were  restricted  on 
assumed  turbulent  flow.  In  Chapter  III,  we  have  developed  a  new  theoretical  framework 
for  prediction  of  the  collision  rate  and  collision  velocity  of  small  particles  at  zero  inertia 
in  general  turbulent  flows.  The  collision  velocity  is  a  critical  variable  related  to  the 
collision  rate.  However,  the  effects  of  particle  inertia  on  the  collision  velocity  and 
collision  angle  have  never  been  predicted  using  Lagrangian  tracking  method.  In  Chapter 
IV,  we  carried  out  an  asymptotic  analysis  to  predict  the  effects  of  small  particle  inertia  on 
the  collision  velocity  and  presented  the  effect  of  the  inertia  on  collision  angle  in  isotropic 
Gaussian  turbulence.  As  pointed  out  in  §1.4,  no  realiable  method  exists  for  quantifying 
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the  concentration  non-uniformity  and  the  effect  of  the  preferential  concentration  on  the 
collision  rate  has  never  been  estimated.  In  Chapter  V,  a  new  quantification  method  has 
been  developed  to  reliably  quantify  the  concentration  non-uniformity  of  particles  in 
turbulent  flows,  and  the  effect  of  the  preferential  concentration  on  the  collision  rate  is 
predicted  in  isotropic  Gaussian  turbulence.  The  first  result  of  Chapter  VI  is  our 
correction  to  the  classical  result  of  Smoluchowski  for  the  collision  rate  of  monodisperse 
particles  in  a  laminar  shear  flow.  The  combined  effects  of  both  the  increased  collision 
velocity  and  increased  concentration  non-uniformity  on  the  collision  rate  due  to  the 
inertia  are  then  examined.  The  effects  of  particle  size,  gravity  and  mean-concentration  on 
the  collision  rate  are  investigated  through  direct  numerical  simulations  in  Chapter  VI. 
The  effects  of  the  collision  models  on  the  collision  rate  are  also  presented  in  Chapter  VI. 


CHAPTER  2 

TURBULENCE  AND  PARTICLE  COLLISION  DETECTION 


2. 1  Turbulence  Representation 
2.1.1  Isotropic.  Gaussian  Turbulence 

The  following  model  for  the  energy  spectrum  developed  in  Mei  &  Adrian  (1995) 
is  assumed  for  turbulence, 

3  k*  1 

B(k) -5  %^[lW],7/e  «„(-nS  k2)  (2.1, 

where  u0  is  the  root-mean-squared  turbulent  velocity,  Icq  is  a  typical  wave  number  and  the 
dimensionless  parameter  Tio^gko  is  related  to  the  turbulent  Reynolds  number  Rex. 
Large  Rex  correspond  to  small  r|0  and  vice  versa.  The  normalizing  coefficient  y  in  E(k) 
is  determined  from  (2.1)  by  satisfying  the  total  energy  requirement, 

l      00        k4  _2 

v"  =    '   =2~^T    exp(-ri0k2)dk  .  (2.2) 

o  (1  +  k  )17/6 

For  small  ti0,  E(k)~k4  when  k=  k/k^l  and  E(k)~k"5/3,  which  is  the  scaling  law  in  the 
inertial  subrange,  when  l«k«l/r|0.  For  large  r\0,  the  above  energy  spectrum  recovers 
that  of  Kraichnan  (1970).  The  relationship  between  the  Eulerian  integral  length  L,,^ 
and  r|0  (or  Rex)  is  given  in  Mei  &  Adrian  (1994). 

A  turbulent  eddy  loses  its  identity  as  it  is  convected  and  dissipated;  this  behavior 
is  described  by  the  eddy-self  decay  function  D(x)  which  has  an  integral  time  scale  TQ. 
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The  Fourier  transformation  of  D(t)  gives  the  power  spectrum  D(co)  .  A  composite  form 
for  the  power  spectrum  D(co)  was  constructed  in  Mei  &  Adrian  (1995)  as 

_  exp(-n.J(o2) 

D(ffl)  =A1^~  (2-3) 

where  r|j  is  intended  for  the  viscous  dissipation  time  scale  if  the  turbulent  Reynolds 
number  is  high,  and  T  is  close  to  T0  if  r|,/T«l.  The  coefficient  A  is  determined  as 

_     t  expC-r^/T2) 

A  =7t  l-erfCrij/T)  (2-4) 

QO  ~ 

by  satisfying  |  D(co)  dto  =1.  The  integral  time  scale  TQ  is  related  to  T  as  TQ  =  7iD(0)  = 


7lA  . 

The  relationship  between  Ln  and  T0  is  not  known  in  general;  it  can  be  taken  to 
have  the  form 

T0  =  c£(Rex)L11/u0,  (2.5) 
where  cE(Rex)  was  determined  approximately  in  Mei  &  Adrian  (1995)  based  on  the 

experimental  data  of  Sato  &  Yamamotto  (1988)  for  the  fluid  dispersion  and  the  value  of 
high  Reynolds  number  turbulent  Prandtl  number. 

Using  random  Fourier  modes  representation,  an  isotropic,  Gaussian,  pseudo 


turbulence 

N 


Uj(x,  t)=  I[b|m)cos(k(m)-x  +co(m)-t  ) 

m=l 

+cjm)  sin(k(m)  •  x  +  co(m)-t  )  ]  (2.6) 
is  constructed  to  simulate  the  turbulent  flow  with  a  specified  energy  spectrum.  In  the 
above  Nk  is  the  number  of  the  random  Fourier  modes  in  one  fluid  realization, 

k(m)and  co(m)  are  the  wavenumber  and  frequency  of  the  m-th  mode.  The  random 
coefficients  b^  and  c^m)  are  chosen  as  follows, 


15 


bf  =bf  (5irk;m)kf  /k(m)  )  f(k(m),  a)(m))  (2.7) 

where  b(m)  follows  a  normal  distribution  with 
J 

<b(^  >  =  0&<b(.m)b(^  >  =  5^  (2.8) 

< 

(m)      im)  (m) 

The  factor  (5{j-kj     kv.    /k     )  in  Eq  (2.26)  ensures  the  incompressibility  Vu=0  for 

every  Fourier  mode.  The  scale  factor  f(k,  co)  depends  on  the  energy  spectrum  and  the 
probability  density  functions  (pdf),  Pi(k)  and  p2(co),  of  k  and  co, 

-2/1      ,  E(k)D(©) 

f  (k,  co)  =  2  — •  (2.9) 

47rN2kPl(k)p2(o)) 

Since  E(k)  decays  as  k~5/3  in  the  inertia  subrange,  a  Gaussian  distribution  for  k  would  lead 
to  a  very  slow  statistical  convergence  for  Uj(x,  t).  To  achieve  convergence  of  the  statistics 

involving  the  derivatives  such  as  — L ,  the  following  algebraic  pdf 

Pii(kj)=         0+  Ikjl)^    for  i=l,  2  &  3  (2.10a) 
and  P,(k)  =  P11(k,)P12(k2)P13(k3)  (2.10b) 

is  used  to  sample  kj  for  1^=1.  The  frequency  co     is  generated  with  the  following  pdf, 

The  scale  factor  is  finally  set  to  be 

A.     ,   E(k)D(co)  

r (k,  (o)  =         2     •  (2. 12) 

4;iNkk  Pn(k1)p12(k2)p13(k3)p2((a) 

It  is  noted  that  in  the  present  studies  on  particle  collisions,  the  collision  rate  is  dictated  by 
the  spatial  structure  and  the  eddy-self  decay  has  no  effect  on  the  collision  rate.  In  the 
direct  numerical  simulation  of  particle  collision  in  isotropic  turbulence,  however,  the 
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eddy  self-decay  can  be  easily  incorporated.  Hence,  the  dependence  of  b .  and  c .  on  co 
is  still  included. 


2. 1 .2  Rapidly  Sheared  Homogeneous  Turbulence 


3\J 

For  an  initially  isotropic  turbulence  under  rapid  shear,  T= — ,  the  variation  of 

dy 

b(.m)  with  time  was  given  in  Townsend  (1976)  based  on  rapid  distortion  theory  (RDT). 

Hunt  &  Carruthers  (1990)  gave  an  overview  on  the  use  of  RDT  for  representing  various 
types  of  homogeneous  turbulence.  Denoting  the  total  shear  S  as 

S  =  Tt  (2.13) 
the  new  wavenumber  vector  x(m)  after  the  shearing  is  given  by 

(m)  _.(m) 

x(m)  =k(m).Sk(jn)  ^ 

(m)  _  .  (m) 
*3  ~K3 

where  k<m)  is  the  wavenumber  vector  of  the  m-th  mode  before  the  mean  shear  is  applied. 
The  amplitude  of  each  Fourier  mode  becomes 

b(m)(S)  =5(m)(())  +aib(™)(0) 

b(2m)(S)  =a2b(2m)(0)  (2.15) 

b(3m)(S)  =b(3m)(0)  +a3b(2m)(0) 

where  b^O)  is  the  corresponding  amplitude  of  the  m-th  mode  before  the  rapid 
distortion.  Dropping  the  superscript  m  for  convenience,  the  coefficient  a{  for  each  mode 
m  is  given  as 

Sk]  k2-2k2VSkik2 
"1=l^k]  X2 
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k2kj        .    ...    k2  ,,  ,k2-Sk, 


(k1+k3)3/2k,  ^k2+k2  ^k2+k2 

a2  =  k2/x2  (2.16) 
Sk,k3  k2-2k2+Sk,k2 


a,  = 


"  k2+k2  X2 


k2k3  .     k2  ,  k2-Skj 


where  k=  ^k2+k2+k2  and  %=  ^Jx^xl+ll  ■  The  coefficient        is  similarly  obtained.  It 

is  clear  from  the  above  that  the  amplitude  of  each  mode  is  entirely  determined  by  the  total 
shear  S=n.  The  dissipation  rate  e(S)  can  be  evaluated  as 

e(S)/v  =  |  <  £  {x<m>2  £  [b(m>2(S)  +cC")2(S)]  }  >.  (2.17) 

m=l  i=l 

Finally,  the  total  velocity  field  is 

Uj(x,t)=  £  [b(.m)(S)  cos(x^-x  +  co(m)t) 
m=l 

+  cf\s)  sin(x(m)-x  +  o)(m)t)] 

+  rx25n.  (2.18) 
Since  the  shear  rate  is  high,  the  term  co(m)t  in  the  above  acts  merely  as  a  random  phase. 


2. 1 .3    Isotropic  Turbulence  Generation  in  a  Periodic  Box 

For  the  direct  numerical  simulations  of  particle  collisions  in  this  work,  the 
velocity  field  of  an  isotropic  and  Gaussian  turbulence  is  generated  by  a  random  Fourier 
modes  representation  (2.6).  To  effectively  conduct  a  computational  study  on  particle 
collision,  a  region  of  finite  size  in  which  a  large  number  of  collisions  occur  should  be 
used  with  periodic  boundary  conditions.  To  render  the  turbulence  represented  by  random 
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Fourier  modes  periodic  in  a  box  of  volume  n3,  the  random  wave  numbers,  k .    ,  are 

rounded  to  the  nearest  even  integers  while  all  other  quantities  are  held  fixed  so  that  the 
flow  is  periodic  within  the  box  of  volume  7t3.  The  periodic  box  of  volume  n3  yields 
virtually  the  same  collision  statistics  when  the  volume  is  increased  to  (2n)3.  Hence  n3  is 
used  throughout  this  study.  The  corresponding  statistics  (turbulence  rms  velocity  uo, 
Reynolds  number  Re^,  Taylor  micro  length  scale  X,  Kolmogorove  length  scale  tj,  and 
Kolmogorov  time  scale  ik)  are  computed  based  on  the  modified  wavenumbers.  As  will 
be  demonstrated  using  the  present  analysis  outlined  in  §2.1.1,  the  collision  statistics, 
when  normalized  using  R3(e/v)1/2  for  an  and  <-wr>,  the  results  are  the  same  using  the 
turbulence  generated  in  an  unbounded  domain  and  in  a  periodic  box. 

Unless  specifically  mentioned,  most  of  the  isotropic  turbulence  used  in  the 
numerical  simulations  throughout  this  study  is  within  the  periodic  box  of  length  n  and 
have  the  following  characteristics:  root-mean-squared  turbulence  velocity  u0  =  1.0, 
integral  length  scale  L, ,  =0.594,  Taylor  micro  length  scale  X=0.217,  Kolmogorov  length 
scale  r|=0.023,  Kolmogorov  time  scale  xk  =  (v/e)l/2  =0.073,  and  Re^=40.1.  The 

turbulence  is  imposed  on  a  box  7i3  with  periodic  boundary  conditions  in  three 
perpendicular  directions.  The  periodic  treatment  of  boundary  conditions  allows  the 
simulation  results  by  a  finite  volume  to  be  extended  to  an  infinite  homogeneous  system. 

Initially,  a  given  number  of  monosize  particles  are  randomly  introduced  in  the 
flow  field  of  volume  n3.  The  periodic  boundary  conditions  are  also  applied  to  the 
motions  of  the  particles.  The  particles  are  advanced  using  Uj(x,  t)  given  by  (2.6)  with  a 
time  step  At.  A  semi-implicit  finite  difference  scheme  is  implemented  to  integrate  the 
equation  of  particle  motion  and  to  detect  collision.  The  scheme  has  an  accuracy  of  second 
order.  The  detailed  error  analysis  of  the  semi-implicit  finite  difference  scheme  will  be 
given  in  next  section. 
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2.2  Numerical  Integration  of  Particle  Motion  and  Error  Analysis 


Major  concerns  in  numerical  simulation  are  the  accuracy  and  robustness  of  the 
numerical  scheme.  Since  the  collision  detection  in  the  present  simulation  is  based  on  the 
calculation  of  particle  trajectories,  the  collision  kernel/rate/velocity/angle  are  affected  by 
the  accuracy  of  the  finite  difference  scheme.  In  this  study,  both  tracking  particles  and 
detecting  collision  are  carried  out  in  a  Lagrangian  frame.  A  semi-implicit  finite 
difference  scheme  is  implemented  to  integrate  the  equation  of  particle  motion  and  to 
detect  collision.  This  scheme  is  unconditionally  stable.  With  consistent  treatment  for 
both  particle  velocity  and  turbulence,  The  truncation  error  in  the  semi-implicit  scheme  is 
0{  (At)2 )  and  independent  of  p.  The  accuracy  analysis  of  the  scheme  is  given  as  follows. 

The  particle  dynamic  equation  in  the  absence  of  (or  weak  )  gravity  is 
dv  n 

dt~=P(u-v)-  (2.19) 
The  semi-implicit  finite  difference  scheme  for  the  particle  dynamic  equation  (2.19)  is 
vn+I-vn  1 

 =  -P[(u-v)n+,+(u-v)n] 

=  |p[(un+,+un-(vn+,+vn)].  (2.20) 
Taylor  series  expansions  for  un+l  and  vn+1  give 


dt     2         dr     6  dt3 


(2.22) 


Each  side  of  equation  (2.20)  can  be  evaluated  as 

dv"     1       d2vn     1  2dV 
-7-  +  T(At)— 2-  +  -(At)2  — 
at     2        dt2     6  dt 


LHS=  — +  -(At)^-  +  i(At)2^-+  O(At)3  (2.23) 


RHS=W«n-yn)4pAt(^-^)  +  Ip(At^(^-^)  +  0(At)3.  (2.24) 
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Then  equation  (2.20)  can  be  rewritten  as 

vn+1-vn  1 

0=  -p[(u-v)n+1  +(u-v)n] 

At  2 

dv"     1       d2vn     1  A  ,2dV 
-(At)— —  +  -(At)2— —  + 


dt     2X  "'  dt2      6V  dt3 

2        dt      dt      4  dt  dt 


Thus,  equation  (2.20)  is  equivalent  to 

-    l^.^.^  .I(At)adV 


2      dt2         dt      dt  ,J    6V    '  dt3 

1        2  d2un  d2vn 
+  -P(At)2(— -^4-). 


Taking  derivative  with  t  repeatedly,  we  get 
and 

Repeatedly  substituting  above  equations  into  the  RHS  of  (2.26)  results  in 
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=  -^-(At)2-^  +  0((At)3).  (2.29) 
12  dr 

It  is  noted  that  the  truncation  error  in  equation  (2.29)  is  independent  of  p.  Thus 
motions  of  particles  with  very  large  P  can  be  simulated  without  loss  of  accuracy. 


2.3  Particle  Collision  Treatment 


2.3.1  Collision  Detection  Scheme 


Collision  detection  between  any  pairing  particles  is  given  as  follows  (Chen  et  al. 
1995).      At  nth  time   step    t0,   these  two  particles   are   initially   located  at 

rio=(xiO'  y  10'  zio)  30(1  r20  =(x20>  v20>  z2o)-  After  a  small  time  advancement  At, 
at  (n+1  )th  time  step  t, ,  new  positions  of  these  two  particles  are  rn  =  (xn,  y , , ,  z, , ) 

and  r2]  =  (x21 ,  y2\ ,  z2] ) .  The  particle  position  vectors  at  time  t,  which  is  between  t0 
and  t  ] ,  can  be  expressed  as 

ri(t)={x10+vxl(t-t0),  y10  +vyl(t-t0),  z10  +v2l(t-t0)}  (2.30) 

r2(t)={x20 +vx2(t-t0),  y20+vy2(t-t0),  z20+vz2(t-t0)}.  (2.31) 
The  distance  between  two  particles,  s(t),  is 

s(t)=|r2(t)-r,(t)|.  (2.32) 
These  two  particles  will  collide  if  the  condition  s(t)<R  is  satisfied.  It  is  noted  that  the 
relative  distance  s(t)  varies  with  time  t,  and  it  is  not  a  monotonic  function  of  t.  The 
minimum  relative  distance  occuring  at  time  tm  is  determined  by 

5s, 


"t  =  t 


m 


=  0-  (2.33) 


Combining  (2.30)-(2.33)  gives 


tm=t0-—  (2.34) 


D2 
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where 

Di  =(x,0  -x20)(vxl  -  vx2)  +  (y10  -y20)(vyl  -vy2) 

+(z10-z20)(vzl-vz2)  (2.35) 

and 

D2  =(vxl  -  vx2)2  +(vyl  -  vy2)2  +(vzl  -  vz2)2.  (2.36) 
Equation  (2.33)  does  not  guarantee  that  the  moment,  when  the  above  mentioned 
minimum  distance  occurs,  has  to  be  within  the  time  interval  t0  and  t0  +  At .  Three  cases 
need  to  be  considered  for  calculating  the  minimum  distance  within  a  time  step. 

(i)  if  tm  <  t0 ,  then  the  minimum  distance  is  at  the  initial  time  t0 . 

smin=s(to)  (2.37) 
i-e-       smin  =  {(x10-x20)2+(y10-y20)2+(z10-z20)2}1/2.  (2.38) 

(ii)  if  tm  >  t0  +  At ,  the  minimum  distance  must  occur  at  the  time  t0  +  At . 

Smin  =s(to+At)  (2.39) 
ie-  smin=(s2+s2+s2)1/2  (2.40) 

where 

sx  =  xio  -x20+At(vxl  -vx2) 
sy  =yio  -yio  +  At(vyl  -Vy2) 

sz  =z10  -Z20  +  At(vzl  ~vz2)- 

(iii)  if  t0  <  tm  <  t0  +  At ,  the  minimum  distance  occurs  during  the  time  interval 
t0  and  t0  +  At .  In  this  case,  we  have 

smin=s(tm)  (2.41) 
ie-  smjn  =  (sL  +  s2^  +s2mz)vl  (2.42) 


where 


Smx  =x10  -x20  +(tm  -t0)(vxl  -vx2) 

Smy  =yiO  ~y20  +(tm  "  t0)(vyl  -  Vy2  ) 
Smz  =z10  "z20  +(tm  -t0)(vzl  -  Vz2). 
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After  smjn  is  evaluated  within  the  time  step,  the  examination  of  collision  events  is 
still  needed  because  smjn  can  be  either  larger  or  smaller  than  R.  Collision  occurs  if 

smin  ^  R  (2.43) 

At  the  collision  instant  tc , 

s(tc)  =  R  (2.44) 
and  s(tc)=|r2(tc)-r1(tc)|.  (2.45) 

Hence  the  collision  instant  tc  can  be  obtained  solving  (2.44)  and  (2.45) 

-b-Jb2-4(s^-R2)w2 

tc  =  t0  +  "  -i9   (2-46) 

2w 

where 

b=2[(x10  -  x20)(vxl  -  vx2)  +  (y10  -y20)(vyl  -  vy2) 

(zio  -z2o)(vzi  -vz2)l 

S0  =  {(XlO  -X20)2  +(yiO  -y20)2  +(ZlO  -z20)2}'/2 

w  =  {(vxl  -  vx2)2  +  (vy,  -  vy2)2  +  (vzl  -  vz2)2}1/2  . 

With  perfect  sticking  collision  model,  any  collision  can  produce  a  new  born  particle 
(daughter  particle).  The  daughter  particle  will  inherit  all  possessions  of  parent  particles, 
i.e.,  conservation  of  mass  and  momentum.  Let  the  daughter  particle  be  numbered  3,  then 
the  conservation  law  gives  the  size  of  daughter  particle  as 

d3=(d?+d2i)1/3  (2.47) 
where  dj  (i=l,  2,  3)  are  particle  radii.  The  velocity  components  of  daughter  particle  are 

vxldl  +  vx2d2 

vx3=  ,3  (2.48) 


d^ 

+  vy2d2 

d^ 

vzld3 

+  vz2d^ 

vy3=  -r^   (2.49) 


vz3=          ,3  (2-50) 

d3 


and  the  position  components  at  the  end  of  time  step 
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x3  =  {[x10  +  vxl(tc  -  t0)]dj  +  [x20  +  vx2(tc  -  t0)]d^}  /  d^  +  vx3(t0  -  tc  +  At)  (2.51) 
y3  =  ([yio  +  vyi(tc  -  to)]d?  +  [y20  +  vy2(tc  -  t0)]d32}  /  d^  +  vy3(t0  -  tc  +  At)  (2.52) 
z3  =  {[zio  +  vzi(tc  "  t0)]di  +  [z20  +  vz2(tc  -  t0)]di}  /  d^  +  vz3(t0  -  tc  +  At)  (2.53) 

In  general,  the  collision  detection  takes  a  lot  of  computer  time  since  a  large 
number  of  particles  are  considered.  The  primitive  collision  detection  scheme  searches  the 
collision  candidates  by  going  through  all  particles  N,  so  that  the  computing  time  is 

proportional  to  N2  .  The  collision  detection  method  in  Chen  et  al.  (1995)  is  an  efficient 
scheme,  and  it  is  also  implemented  in  this  thesis.  By  dividing  the  computational  domain 
into  a  number  of  small  cells,  the  potential  collision  partners  for  a  given  particle  in  one  cell 
are  then  searched  within  this  cell  and  neighboring  twenty-six  cells  during  one  time  step. 
This  search  only  involves  a  small  number  of  particles.  The  total  computing  time  with  this 
scheme  can  be  reduced  dramatically.  Only  binary  collision  is  considered  by  assuming  a 
negligible  probability  of  multiple  collisions  in  a  small  time  step  in  a  dilute  condition.  For 
any  particle  moving  out  of  the  computational  box,  it  is  re-introduced  into  the  box  by 
invoking  periodicity.  The  collision  rate  or  kernel  can  be  determined,  in  principle,  by 
counting  the  number  of  collisions  per  time  step.  Initially,  a  given  number  of  particles  are 
randomly  distributed  into  the  computational  domain.  The  particles  are  evolved  with 
turbulence  for  a  period  of  time  before  collision  detection  turned  on.  The  collision  events 
are  not  counted  until  all  particles  with  finite  inertia  have  reached  dynamical  equilibrium 
in  turbulence. 

2.3.2  Post-Collision  Treatment 

The  collision  event  is  counted  when  two  particles  are  brought  into  contact,  as 
stated  above.  The  collision  among  particles  and  its  subsequent  treatment  are  expected  to 
disturb  the  original  system  and  permanently  impact  the  subconsequent  collision 
evaluation.  For  the  droplets  in  gas  turbulence,  the  collision  may,  most  likely,  result  in 
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coalesce  of  two  droplets  and  a  larger  droplet  is  then  formed  by  the  colliding  droplets. 
Similar  phenomena  may  be  found  in  a  system  consisting  of  bubbles  and  fluids.  The 
collision  between  two  solid  particles  may  behave  in  different  ways.  Two  solid  particles 
may  stick  together  upon  contact  or  the  colliding  solid  particles  will  separate  after 
collision.  In  summary,  the  issue  dealing  with  immediate  consequence  of  those  contacts  is 
referred  as  post-collision  treatment.  Several  models  with  regard  to  post-collision 
treatment  have  been  developed  during  past  decades. 

In  Chen  et  al.  (1995)  every  collision  results  in  the  collision  partners  disappearing 
from  their  respective  size  groups  to  produce  a  larger  particle.  Mass  and  momentum 
conservation  laws  were  applied  for  the  birth  of  the  new  particle.  This  scheme  is 
considered  to  be  a  physically  realistic  model  assuming  a  unity  collision  efficiency.  It  is 
hereinafter  referred  as  'KEEP  ALL'  model.  Balachandar  (1988)  used  a  different  post- 
collision  scheme.  Starting  with  monosize  particles  in  system,  after  each  collision,  the 
resulting  larger  particle  was  discarded  in  the  computation  so  that  it  does  not  contribute  to 
future  collisions.  This  post-collision  treatment  is  hereinafter  referred  as  the  'THROW 
AWAY'  scheme.  Sundaram  &  Collins  (1997)  employed  a  'HARD  SPHERE'  collision 
model  in  which  the  particles  were  assumed  to  be  rigid  spheres  and  forced  to  bounce  back 
after  the  contact.  Momentum  conservation  is  applied  to  the  particles  after  collision.  In 
contrast  to  the  'KEEP  ALL'  and  'THROW  AWAY'  schemes,  the  'HARD  SPHERE' 
model  allows  a  particle  to  have  multiple  collisions  with  other  particles  in  a  system.  Wang 
et  al.  (1997)  tested  three  different  schemes  and  presented  one  of  their  own.  In  their 
scheme,  collision  does  not  disturb  the  system.  The  colliding  particles  are  allowed  to 
overlap  in  space  (hereinafter  referred  as  'OVERLAPPING'  scheme),  but  no  larger 
particles  are  formed  after  the  collision.  Two  colliding  particles  can  separate  after  their 
trajectories  satisfy  |r2(t)  -  r,(t)|>R.  Thus,  a  particle  can  have  multiple  collisions  with 
other  particles.  These  four  different  post-collision  treatments  will,  in  principle,  result  in 
different  collision  rate/kernel. 
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Figure  2. 1    The  effect  of  time  step  on  the  cumulative  average  of  the  normalized 
collision  kernel  an  for  inertialess  particles  in  turbulence.  Particle  size  D=0.05. 


2.4  Computer  Simulation 

2.4.1  Time  Step 

Three  collision  models:  'THROW  AWAY,"  "KEEP  ALL"  and  "HARD 
SPHERE"  are  employed  in  this  work.  For  the  following  simulations  with  "THROW 
AWAY"  and  "KEEP  ALL"  collision  models,  the  time  step  At=0.01  is  used  to  advance  the 

motions  of  particles.  Figure  2.1  shows  the  effect  of  time  step  on  the  collision  kernel  an 
using  "KEEP  ALL"  collision  model.  It  can  be  seen  that  there  is  almost  no  difference 
between  the  results  of  At=0.01  and  At=0.005.  Similar  results  are  obtained  by  'THROW 
AWAY"  scheme.  The  volume  concentration  of  suspended  particles  is  usually  quite  low 
in  many  natural  systems  such  as  spray  burners  and  air  scrubbers,  which  is  within  the 
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interest  of  this  work.  In  such  dilute  systems,  the  multiple  collisions  during  a  small  time 
step  may  be  negligible,  thus  binary  collision  mechanism  is  considered  solely  in  this 
study. 

2.4.2  Turbulence  Realization 

The  averages  of  statistical  quantities  are  evaluated  in  two  directions:  average  over 
all  time  steps  and  turbulence  ensemble  average  over  all  turbulence  realizations.  The  total 
time  steps  mainly  depend  on  the  particle  response  time  and  the  number  of  collision  pairs 
during  whole  time  period.  We  start  counting  the  time  steps  for  collision  detection  after 
the  motions  of  particles  have  been  evolved  for  three  times  the  particle  response  time. 
Ideally,  to  form  the  turbulence  ensemble  average,  the  numerical  simulation  should  be 
repeated  an  infinite  number  of  realizations  with  different  initial  particle  distributions.  An 
infinite  number  of  turbulence  realizations  is  impractical  due  to  limited  computer  capacity 
and  may  be  unnecessary  from  the  viewpoint  of  the  accuracy  of  numerical  scheme. 
Figures  2.2  (a)  &  (b)  show  the  turbulence  dissipation  rate  against  turbulence  realizations. 
While  the  ensemble  averaged  turbulence  dissipation  rate  varies  with  the  realization  by  a 
large  amplitude  shown  in  figure  2.2  (a),  the  cumulative  average  of  turbulence  dissipation 
rate  approaches  a  constant  for  the  number  of  turbulence  realizations  greater  than  30,  as 
seen  from  figure  2.2  (b).  Figures  2.3  (a)  &  (b)  show  the  small  eddy  shear  rate  against 
turbulence  realizations.  It  can  also  be  seen  that,  although  the  ensemble  averaged  small 
eddy  shear  rate  randomly  varies  with  the  turbulence  realization  in  figure  2.3  (a),  the 
cumulative  average  of  small  eddy  shear  rate  is  close  to  a  constant  by  the  40th  realization 
seen  in  figure  2.3  (b).  Therefore,  in  this  study  forty  turbulence  realizations  are  used  for 
turbulence  ensemble  average. 

2.4.3  Validation  of  the  Collision  Detection  Scheme 

Starting  with  monosize  particles  in  system  and  using  the  'THROW  AWAY" 
collision  model  for  post-collision  treatment.  The  evolution  of  the  particle  number  n^t) 
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follows  the  simple  population  balance  equation: 

dn^dt  = -annf       with     n,(t  =  0)  =  n0.  (2.54) 
The  solution  is 

n0/n,(t)  =  l  +  a„n0t.  (2.55) 
It  should  be  mentioned  that  the  collision  rate  of  monosize  particles  (2.54)  is  only  accurate 
for  large  number  of  particles  in  a  system.  For  small  number  of  particles,  however,  the 
collision  rate  (2.54)  must  be  corrected.  The  detailed  discussions  on  this  issue  will  be 
given  in  §6.2.  Since  the  statistics  for  the  number  of  the  remaining  particles  n^t)  is  much 
more  accurate,  the  average  slope  of  the  curve  1  /  rii(t)  gives  the  collision  kernel,  so  it  can 
be  used  to  check  the  accuracy  of  the  collision  kernel  obtained  by  simulation  in  laminar 
shear  flow. 

To  validate  the  collision  detection,  particle  collision  in  a  uniform  shear  flow  of 
velocity  gradient  T=1.0  is  first  considered.  For  large  number  of  particles  the  theoretical 
prediction  of  the  collision  kernel  ay  of  Smoluchowski  is  accurate  and  is  given  by 

a    =  4  R  ?;  T  /  3  .    For  a  small  number  of  monosized  particles  the  theoretical  prediction 

ij  u 

of  the  collision  rate  of  Smoluchowski  has  to  be  corrected,  as  mentioned  above.  Ten 
thousand  particles  of  radius  0.075  are  introduced  into  a  box  of  volume  2x2x2.  The  time 
step  is  Dt=0.01  and  a  total  of  800  time  steps  is  used.  The  number  of  shear  flow 
realizations  is  80.  Since  no  periodic  condition  may  be  imposed  along  vertical  direction  in 
uniform  shear  flow,  periodic  boundary  conditions  are  employed  in  the  streamwise  and 
spanwise  directions  only.  In  such  a  case,  a  boundary  correction  of  the  collision  rate  has 
to  be  imposed  (Wang  et  al.  1997).  Another  correction  on  the  collision  rate  is  also  applied 
due  to  the  inaccurate  expression  of  the  collision  rate  (§6.2)  Figure  2.4  shows  the 
particle  number  evolution.  The  collision  kernel  a  [  ]  can  be  obtained  from  the  average 

slope  of  the  curve.  Very  good  agreement  can  be  observed  between  simulation  result  and 
Smoluchowski 's  prediction.  Hence,  the  collision  detection  scheme  and  the  numerical 
implementation  are  validated. 
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Figure  2.2  Turbulence  dissipation  rate  s  against  turbulence  realization,  (a)  Variation  of 
ensemble  average,  (b)  Cumulative  average. 
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Figure  2.3  Turbulence  small  eddy  shear  rate  (e  /  v)  against  turbulence  realization, 
(a)  Variation  of  ensemble  average,  (b)  Cumulative  average. 
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Figure  2.4  Inverse  of  total  number  of  monodisperse  particles  against  time  in  laminar 
shear  flow.  T=1.0,  D=0.015. 

For  "HARD  SPHERE"  particles,  the  time  step  for  advancing  the  motions  of 
particles  is  much  less  than  0.01.  It  is  found  that  At=0.001  may  generate  almost  same 
result  as  At=0.0005  for  large  inertia  particles.  For  small  inertia  or  zero  inertia  particles, 
the  collision  kernel  strongly  depends  on  time  step  due  to  multiple  collision  mechanism, 
which  is  shown  in  figures  2.5  for  uniform  shear  flow.  A  pair  of  particles  can  have  tens  or 
hundreds  of  collisions  depending  on  the  size  of  time  step  before  they  completely  separate. 
Similar  multiple  collisions  have  been  found  for  particles  in  turbulence.  Due  to  this 
technical  difficulty,  the  "HARD  SPHERE"  collision  model  is  limited  to  particles  with 
large  inertia. 
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Figure  2.5  Multiple  collisions  of  HARD  SPHERE  particles  in  laminar  shear  flow. 


CHAPTER  3 

A  NEW  THEORETICAL  FRAMEWORK  FOR  PREDICTING  COLLISION 
RATE  OF  SMALL  PARTICLES  IN  GENERAL  TURBULENT  FLOWS 


A  new  theoretical  framework  is  developed  to  predict  the  collision  rate  and  collision 
velocity  of  small,  inertialess  particles  in  general  turbulent  flows.  The  present  approach 
evaluates  the  collision  rate  for  small  particles  in  a  given  instantaneous  flow  field  based  on 
the  local  eigen-values  of  the  rate  of  strain  tensor.  An  ensemble  average  is  applied  to  the 
instantaneous  collision  rate  to  obtain  average  collision  rate.  The  collision  kernels 
predicted  by  Smoluchowski  (1917)  for  laminar  shear  flow  and  by  Saffman  &  Turner 
(1956)  for  isotropic  turbulence  are  recovered.  The  collision  velocities  presently  predicted 
in  both  laminar  shear  flow  and  isotropic  turbulence  agree  well  with  the  results  form  direct 
numerical  simulation  for  particle  collision  in  both  flows.  The  present  theory  for 
evaluating  the  collision  rate  and  collision  velocity  is  also  applied  to  a  rapidly  sheared 
homogeneous  turbulence.  Using  the  mean  turbulent  shear  rate  (s/v)  as  the 
characteristic  shear  rate  to  normalize  the  collision  rate,  the  effect  of  the  turbulence 
structure  on  the  collision  rate  and  velocity  can  be  reasonably  described.  The  combined 
effects  of  the  mean  flow  shear  and  the  turbulence  shear  on  the  collision  rate  and  collision 
velocity  are  elucidated. 


3.1  Introduction 

Smoluchowski  (1917)  considered  collision  rate  among  spherical  particles  in  a 
laminar  shear  flow,  (ux,  ly  uz)  =  (fy,  0,  0)  with  a  constant  gradient  Y.   For  a  target 
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particle  of  radius  rx  centered  at  an  arbitrary  position  xi?  any  particle  with  a  radius  rj 
moving  toward  the  target  particle  will  cause  a  collision  if 

|Xi  -  xjl  <  Kg  =  vt  +rj  (3.1) 
is  satisfied.  The  collision  rate  can  be  evaluated  as 


nc 


..  =  n;nj  J    ft-wrdA  (3.2) 

ij     ^  J    wr<0  r 

where  wr  is  the  radial  component  of  the  relative  fluid  velocity  between  X;  and  Xj,  wr<0 
indicates  an  impending  collision,  -n,n2wr  is  the  particle  flux  moving  inward  to  the  target 
particle,  and  A  is  a  spherical  surface  with  radius  Rjj=rj  +r:.  Using  standard  spherical 
coordinates  (R,  9,  <)>)  centered  at  particle  i  (see  figure  3.1),  wr  on  the  spherical  surface  of 

2 

R=Rjj  can  be  written  as  wr  =  TycosG  =  Ry  sin0  cosG  cos(j).  With  dA=  Ry  sinG  dGd(j),  it  is 


easily  seen  that 


4  3 

nCjj  =  njnj  3  T  Ry  (3.3) 


or 


Oy-frRj.  (3.4) 

Saffman  &  Turner  (1956,  hereinafter  referred  as  ST)  presented  a  classical  theory 
for  the  collision  rate  of  small  droplets  in  Gaussian,  isotropic  turbulence.  The  collision 
rate  among  particle  size  group  i  and  j  is 

nCij  =  riiHj  <  JWr<Q  -wrdA  (3.5) 

The  ensemble  average,  denoted  by  <  >,  is  necessary  since  turbulence  is  random.  Because 
of  the  continuity  of  the  fluid  flow,  the  volume  influx  entering  the  spherical  surface  of 
radius  R^  is  equal  to  the  efflux  so  that 

nc.=^njni<J   .     .     |wr|dA>.  (3.6) 

"c,j     2    ^  J       entire  sphere  1    r'  v  ' 

ST  further  interchanged  the  integration  with  the  ensemble  average  so  that  nCy  can  De 
evaluated  in  isotropic  turbulence, 
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nc..*i  nm:  \        .     <|wl>dA.  (3.7) 

CU     2    ^  J    entire  sphere    1    r'  v  ' 


Figure  3.1  Local  coordinates  for  collision  rate  calculation  in  laminar  shear  flow. 


For  isotropic  turbulence,  the  above  becomes 

nCjj  =  n^  27tRy  <|wx|>  (3.8) 

where  x  is  the  local  coordinate  parallel  to  the  line  connecting  the  centers  of  the  colliding 
particles  as  shown  in  figure  3.2.  For  particle  size  smaller  than  Kolmogorove  length  scale 
T|,  the  relative  velocity  upon  the  collision  was  evaluated  using  the  standard  results  of 
small-scale  turbulence  theory  as 

<KI>  =  <|uix-  ujx|>  *  Ru<^>  =  Rij(]|^^)1/2  (3-9) 

where  e  is  the  dissipation  rate  of  turbulence  and  v  is  fluid  kinematic  viscosity.  This  gives 

n^n^  [yj -] 1/2  =1.29441^ Q 1/2  (3.10) 

or  aij=1.2944R?j(^)1/2  (3.11) 

or  aj  =— 71—  =1-2944  (3.12) 
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where  a::  is  normalized  collision  kernel. 


Figure  3.2  Sketch  of  particle  collision  velocity.  The  separation 
distance  is  X2-X1  =  ±  Ry=  ±(rj+rj). 

It  needs  pointed  out  that  from  equation  (3.6)  to  equation  (3.7),  the  interchange 
between  the  integration  with  the  ensemble  average  is  only  an  approximation  even  in  a 
Gaussian,  isotropic  turbulence.  Furthermore,  the  interchange  between  these  two 
operations  severely  limits  the  extension  of  ST's  result  and  approach  to  other  turbulent 
flows  such  as  the  near  wall  region  where  the  mean  flow  gradient  contributes  to  <|w  |>.  Of 

course,  in  the  theory  of  Smoluchowski  for  laminar  shear  flow,  this  issue  does  not  arise 
since  the  ensemble  average  is  not  needed.  In  a  typical  industrial  facility  where  particle 
collision  leads  to  the  desirable  growth  of  particles,  the  inhomogeneity  of  the  turbulent 
flow  structure  is  common.  Due  to  the  effects  such  as  turbopheresis,  particles  tend  to  be 
driven  toward  wall  regions  to  cause  a  further  inhomogeneity  in  the  particle  concentration 
distribution.  It  was  not  clear  how  one  can  predict  the  collision  rate  of  small  particles  in 
the  presence  of  both  strong  mean  flow  shear  and  high  turbulence  intensity  (Kontamaris, 
1995). 

In  this  chapter,  a  new  theoretical  framework  is  presented  for  the  collision  rate  of 
inertialess  particles,  whose  sizes  are  smaller  than  the  Kolmogorove  length  scale,  in 
general  turbulent  flows.  The  present  analysis  starts  from  equation  (3.5).  Using  the 
leading  order  term,  (xj-xj)-Vu,  in  the  Taylor  series  expansion  to  approximate  the  relative 
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fluid  velocity,  w,  at  the  collision  instant,  the  collision  rate  for  an  arbitrarily  given  Vu  is 

evaluated  first.    The  ensemble  average  is  then  carried  out  in  a  Gaussian,  isotropic 

turbulence  to  obtain  the  collision  rate  and  collision  velocity.     Direct  numerical 

simulations  are  also  performed  to  obtain  the  collision  rate  and  collision  velocity  in  a 

Gaussian,  isotropic  turbulence.    Good  agreement  is  obtained  between  the  present 

* 

prediction  and  the  simulation  for  the  collision  kernel  .  Excellent  agreement  is 
obtained  for  the  average  collision  velocity  between  the  prediction  and  the  simulation. 
While  the  present  result  for  the  collision  rate  agrees  with  the  prediction  of  ST,  the  present 
theory  predicts  an  average  collision  velocity  <-w,>  that  is  1.58  times  the  value  given  by 
ST.  The  difference  is  caused  by  the  bias  of  the  collision  toward  higher  values  of  collision 
velocity  which  was  not  considered  in  previous  studies  on  flow  induced  particle  collisions. 
Using  the  rapid  distortion  theory  (RDT)  for  a  rapidly  sheared  homogeneous  turbulence 
and  the  random  Fourier  modes  representation,  the  dependence  of  particle  collision  kernel 
ay  on  the  mean  flow  shear  rate  and  the  turbulence  structure  in  homogeneous  shear  flow  is 
studied  under  the  present  theoretical  framework.  It  is  found  that  the  collision  kernel,  after 
normalized  by  particle  volume  and  the  turbulence  mean  shear  rate  characterized  by  the 
dissipation  rate,  mainly  depends  on  the  ratio  of  the  mean  flow  shear  rate  to  the  turbulence 
mean  shear  rate. 


3.2  Analysis 

The  relative  fluid  velocity,  w  =  u,  -u2,  between  target  particle  centered  at  X!  and 
the  colliding  particle  centered  at  x2  in  space  with  |x,  -  x2|  =  |r12|  =  R  less  than  the 
Kolmogorove  length  scale  r\  can  be  expressed  using  Taylor  series  expansion  as 

w  =  u2  -u,  *r12-Vu  +  ...  (3.13) 
r12=   x2-x,  (3.14) 

to  the  leading  order.  The  collision  rate,  nq2  ,  for  a  given  Vu  at  a  given  instant  is  thus 
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nq2  =  nin2  lWr<0  -w-ndA  (3.15) 

where  n  is  the  outward  normal  of  the  surface  given  by  |x2  -  X\\  =  Ry  =  R  and  wr=  w  n. 
For  a  spherical  surface,  n=  r12/R  so  that 


\  _  -w-ndA  =  o  J 


I  -WIlU/\  —  D   J       „  -Wi  n 

wr<0  R    w„<0  u 


•w-r,-)  dA 


*R'wr<0-r>2,V,,'r12dA-  (3-16) 

Representing  Vu  as  the  sum  of  a  symmetrical  part  (rate-of-strain  tensor)  and  an  anti- 
symmetrical  part  (vorticity  tensor),  it  can  be  easily  seen  that  the  antisymmetrical  part  of 
Vu  does  not  contribute  to  wr.  Physically,  it  is  simply  because  a  rigid  body  rotation  does 

1  T 

not  contribute  to  normal  velocity.  For  a  symmetric  tensor  E  =  2  +(^u)  )  >  or  a 
symmetric  matrix  Ey,  a  linear  transformation,  which  is  a  rotation  of  the  coordinate  system 
from  (x,  y,  z)  to  (x  ,  y  ,  z  ),  can  be  easily  found  to  reduce  Ey  to  a  diagonal  form  as 


E  = 


ex  0  0 
0  ey  0 


(3-17) 


0  0  e2 

where  (ex  ,  e^  ,  e2  )  are  the  eigen-values  of  E  and  are  the  principal  values  of  the  rate-of- 
strain  tensor.  Due  to  incompressibility,  ex  +  ey  +  e2  =  0  so  that  there  are  only  two 

independent  parameters  (say  ex  and  ey  )  in  E  in  order  to  evaluate  the  integral  in  (3.16). 
Arranging  the  eigenvalues  in  the  order  of 

e2  =-(e*  +ey)%  <ex,  (3.18) 

and  defining 

ey 

C,  =  T-  with    -0.5<£<1  (3.19) 

the  integrand  in  (3.16)  can  be  expressed  as 

-  rI2  -Vu-i-,2  =  -  (ex  x  2  +  e?  y  2  +  e2  2  2) 

—  «k  [x2  +  ;y2-d+Q  22] 
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=  -ex  R2[cos20+;sin20cos2f(l^)sin20sin2<(>]  (3.20) 

where  a  local  spherical  coordinate  system  (x  =Rcos0,  y  =Rsin0cos<|),  z  =Rsin0sin<(>)  has 
been  introduced.  Hence 

r  *  2nTC 

Jw  <Q  -wndA=  -R3ex  ffw  <Q[cos20+  £sin20cos2<t>-(l+Osin20sin2(|>]  sin0  d0  d<|> 
0  0 

=  R3exfl(Q.  (3.21) 
Again,  wr<0  in  the  integration  limit  implies  that  only  the  non-positive  values  inside  the 
square  bracket  are  evaluated  in  the  integration.  While  an  analytical  expression  for  HQ  is 
difficult  to  obtain  due  to  the  dependence  of  the  integration  limits  on  arbitrary  £,  accurate 
numerical  integration  for  HQ  can  be  easily  carried  out  (using  Simpson's  rule).  The 
following  interpolation  is  then  constructed  for  HQ  based  on  numerical  integration 
results, 

HQ  *  f  [0.90725  +  0.2875(0.5+Q2  +  0.333(0.5+Q4]  £<0 

*|  [1  +  0.55£  +  0.72407 1^2  -  0.96874£3 

+  0.73924C4  -  0.23035;5].  C>0  (3.22) 

r  8 

Consider  the  laminar  shear  flow  for  example.  In  this  case,  Q=0,  ex  =  y ,  #(0)  =  j ,  so  that 
18  4 

nc,2  =  nin2  2       3  =  nin2  3       w°ich  is  the  same  as  given  by  Smoluchowski  (1917). 

Equation  (3.22)  is  valid  for  -l<£<-0.5  as  well  since  it  is  based  on  (3.21).  The  restriction 
for  £>-0.5  results  only  from  equation  (3.18)  when  the  eigenvalues  are  arranged  in  such  a 
manner.  For  -1<<^<0,  HQ  is  symmetrical  with  respect  to  C,  =-0.5.  For  example,  £=0 

corresponds  to  a  2-D  stagnation  flow  in  (x  ,  z  )  plane  and  Q=-\  corresponds  to  a  2-D 

stagnation  flow  in  (x  ,  y  )  plane  so  that  the  collision  rates  are  completely  equivalent.  It  is 
also  interesting  to  note  that  £=-0.5  is  equivalent  to  £=1  since  both  correspond  to  an 
axisymmetric  stagnation  flow  with  £=-0.5  for  an  axisymmetric  contraction  and  Q=\  for  an 
axisymmetric  expansion;  the  only  difference  is  the  scale  ex  .   With  these  examples  in 
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mind,  a  simple  analysis  for  the  flow  field  quickly  shows  that  #{£<0)  can  be  obtained  from 
!H{C>0)  as 

Hr^  )  =  -1)  =       HQ         for      <;>0.  (3.23) 

Finally,  the  collision  rate  for  small  inertialess  particles  in  a  turbulent  flow  is  simply 

ncl2  =  n1n2R3<ex  fl(Q>  (3.24) 

where  the  ensemble  average  is  taken  over  all  possible  values  of  the  rate-of-strain  tensor. 
It  is  emphasized  that  this  result  is  applicable  to  general  turbulent  flow.  Needless  to  say, 
the  isotropy  assumption  is  not  needed  in  the  present  framework.  When  the  detailed 
knowledge  of  Vu  is  available,  such  as  from  DNS  or  random  Fourier  modes 
representation  (Kraichnan,  1970)  of  turbulence,  the  ensemble  average  in  (3.24)  can  be 
carried  out. 

The  foregoing  analysis  can  also  be  extended  to  evaluate  the  average  of  the  normal 
(or  radial)  component  of  the  collision  velocity  <-Wf>.  For  a  given  position  with  arbitrary 
ex  and  C,  at  any  given  instant,  the  instantaneous  average  of  the  normal  collision  velocity 
can  be  calculated  by  integrating  over  the  spherical  surface  of  radius  R  as 

<wr|ex,^=  f  .  (3.25) 

wr<0  r 

In  the  above,  the  superscript  Jl  on  the  left  hand  side  denotes  the  area-average,  (-wr)dA  in 
the  numerator  can  be  interpreted  as  proportional  to  the  probability  of  a  particle  striking 
the  target  surface  at  a  position  (0,  <)>)  over  an  area  dA  and  wr  in  the  numerator  is  the 
relative  radial  velocity  of  the  collision  at  (0,  <(>)  on  the  target  surface.  The  denominator  is 
the  normalizing  factor  for  the  probability.  The  product,  wr(-wr),  clearly  indicates  that 

<wrlex  >  is  a  biased  average  in  favor  of  these  regions  in  (0,  <)))  where  (-wr)  is  high. 
Suppose  that  NR  realizations  of  turbulence  be  taken  in  a  turbulent  flow,  the  total  number 
of  collisions  will  be 

Nc  =  NRn,n2</Wr<0  (-wr)dA> 
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per  unit  time  using  equation  (3.5).  The  chance  of  a  particle  striking  the  target  surface 

over  an  area  dA  at  (0,  <|>)  relative  to  the  axis  of  the  principal  direction  x  of  the 
instantaneous  rate-of-strain  tensor  in  a  given  realization  is  then  nin2(-wr)dA|  <Q/NC. 

Integrating  over  the  target  surface  and  summing  over  NR  realizations,  the  ensemble 
average,  <w,>,  is  obtained, 

f     ^  <f  „(-wr)w,dA> 

nin2(-wr)  wr<ov  ry 


<  Wj>  =  Nr  <  f     wr-^  dA>  =    <f'     (    ^  .  (3.26) 

wr<0 

The  above  can  be  simplified  as 


where 

#Q=2ff    [cos2e+;sin29cos2<t)-(l+^)sin2esin2<))]2       sinG  d0  dty.  (3.28) 
0  0 


The  integration  can  again  be  evaluated  numerically  and  represented  in  the  following 
piecewise  interpolations 

^  =0.4+1. 0777(0.5+Q2  -  0.64767(0.5+Q4         for  <;<0 

=0.62831  +  0.68495^  -  0.09 1083^2 

+  0.09030 1<;3  -  0.034658^4  for      <;>0.  (3.29) 

Similarly,  the  mean  square  value  of  the  radial  component  velocity  <  wr  >  can  be 
evaluated  as 

2  <e*«GK)> 
<W'  >  =  R<e^Q>  (3-3°) 

with 

^  =0.171 14+1.2643(0.5+Q2-0.48123(0.5+Q4    for  <;<0 
=0.45714  +  0.97655<;  +  0.39606C2 
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+  0.072430C3  -  0.029 1£4  for      C>0.  (3.31) 

The  functions  #[Q,  §(Q/fl(Q,  and  flJQ/^Q  are  shown  in  figure  3.3.  Take  a  laminar 
shear  flow  for  example:  ex  =172  and  £=0  so  that  0.62832  (=rc/5)  and  £/#=16/35  . 

It  is  easily  seen  that 

„  r2  7i  8  r  8 ,     7i  _ 

<wr>  =  -R^53/[23]  =  -To  TR,  (3.32) 

2      ,  r3 16  8  r  8     4  , 

<Wr  >  =  R  ¥35  3  ^23^35™   •  (333> 

The  results  of  these  two  moments  of  wr  can  be  used  to  validate  the  present  theoretical 
formulation  by  comparing  with  direct  numerical  simulation  for  particle  collisions. 


-0.5  0.0  0.5  1.0 


Figure  3.3      Variations  of  HQ,  #Q,  KLQ  with  S  where 
f.-8  ^Q,f2-^gandf3=  — . 
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To  apply  the  above  results  to  turbulent  flows,  we  consider  the  following  two 
cases:  i)  Gaussian  isotropic  turbulence;  and  ii)  rapidly  sheared  homogeneous  turbulence 
under  a  mean  flow  gradient  T.  In  both  cases,  the  turbulence  can  be  represented  using 
random  Fourier  modes  and  ensemble  averages  can  be  obtained  without  resorting  to  the 
direct  numerical  solutions  to  the  Navier  Stokes  equations  for  the  fluid  flows.  Pertinent 
comparison  can  be  made  with  ST's  prediction  in  the  first  case.  Insights  on  the  effect  of 
mean  flow  shear  can  be  gained  by  examining  the  collision  rate  in  the  rapidly  sheared 
turbulence. 

3.3     Numerical  Simulation 

The  detailed  representation  of  the  isotropic  turbulence  and  rapidly  sheared 
homogeneous  turbulence  are  given  in  Chapter  2.  In  short,  the  isotropic  turbulence  is 
represented  in  the  form  of 

Ui(x,  t)=  I[bSm)cos(k(m) -x  +co(m)-t  ) 

m=l 

+cSm)sin(k(m)-x  +w(m)-t  )  ]  (2.6) 

where  k(m)and  co(m)  are  the  wavenumber  and  frequency  of  the  m-th  mode.  The  random 

coefficients  b-m)  and  c-m)  are  their  corresponding  amplitudes.  For  sheared 
homogeneous  turbulence,  the  flow  field  is  given  by 

Uj(x,  0=  Z[b|m)(S)cos(X(m)-x+co(m)-t  ) 

m=l 

+  cSm)(S)sin(x(m)  -x+co(m)  -t  ) 

+  rx25i,  (2.18) 
where  x<m)  is  the  new  wavenumber  after  the  distortion  by  the  mean  flow  shear,  S  =  Tt  is 
the  total  shear  by  the  mean  flow,  and  b^VS)  is  the  amplitude  of  the  m-th  random  mode 
after  the  rapid  shearing. 
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The  scheme  developed  in  Chen  et  al.  (1995)  is  employed  here  to  search  for 
particle  collisions.  Detailed  description  on  the  collision  algorithm  is  given  in  Chapter  2. 
Periodic  boundary  conditions  are  employed  for  both  turbulent  flow  and  particle  motion.  A 
number  of  particles  (Np)  are  randomly  introduced  in  the  flow  field  of  volume  tc3.  The 
average  collision  velocity  is  obtained  based  on  all  the  collision  events  for  whole  time 
period  and  all  turbulence  realizations. 

As  discussed  in  Chapter  2,  four  different  post-collision  treatments  will  give 
different  collision  rate  and  it  is  not  clear  which  treatment  is  the  best.  In  the  study 
throughut  this  chapter,  the  'THROW  AWAY'  scheme  is  implemented.  However,  an 
extremely  low  particle  volume  concentration  is  used.  The  initial  particle  volume 
concentration  is  chosen  carefully  so  that  the  average  number  of  collisions  in  each 
turbulence  realization  is  less  than  1  .  The  post-collision  treatment  thus  has  little  practical 
effect  on  the  collision  statistics.  A  large  number  of  turbulence  realizations  are  used  to 
obtain  sufficient  ensemble  average.  A  subsequent  cumulative  time  average  smoothes  out 
statistical  noise  in  the  collision  rate. 


3.4  Results  and  Discussions 
3.4. 1    Comparison  with  Simulation  Results  in  a  Laminar  Shear  Flow 

In  an  earlier  paper  dealing  with  particle  collision  a  laminar  shear  flow  (Hu  &  Mei, 

4  ^ 

1998),  comparison  for  the  collision  kernel,  j  TR  ,  based  on  Smoluchowski's  prediction 

and  the  direct  numerical  simulation  was  presented;  excellent  agreement  was  obtained.  To 
confirm  the  present  analyses  in  the  laminar  shear  flow,  we  compare  the  collision  velocity 
given  by  equations  (3.32-3.33)  with  that  based  on  the  direct  numerical  simulations.  Zero 
inertial  particles  of  radius  0.009,  0.0125  and  0.02,  respectively,  are  introduced  into  a 
cubic  box  with  length  Lx=2,  Ly=2.0,  and  Lx=2.0,  with  the  restriction  0.1<yp<1.9.  A 
uniform  shear  flow  with  a  shear  rate  T=l  .0  is  imposed  and  the  flow  velocity  field  is  given 
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by  v=(u,  v,  w)  =  (ry,  0,  0).  Periodic  boundary  conditions  are  employed  in  the  streamwise 
and  spanwise  directions.  Particles  near  y=0.1  and  y=1.9  have  less  collisions  because 
there  are  no  particles  in  regions  of  y<0.1  and  y>1.9  in  the  simulation.  A  correction  due  to 
this  boundary  effect  (Wang  et  al.  1997)  was  needed  for  the  collision  rate  (Hu  &  Mei, 
1998).  However,  this  boundary  correction  has  little  effect  on  the  collision  velocity  since 
the  decrease  in  the  collision  kernel  affects  both  the  numerator  and  denominator  in 
equation  (3.25).  Three  hundred  time  steps  with  a  step  size  At=0.01  are  used  to  advance 
particles  and  a  total  of  40  realizations  were  employed  to  further  increase  the  statistical 
accuracy.  Ten  thousand  particles  are  used  for  r=0.009,  0.0125  while  four  thousand 
particles  are  used  for  larger  size  with  r=0.02.  Table  3.1  shows  the  comparisons  for  both 

<-w,>  and  <wf  >.  Since  <wf  >  is  the  second  order  moment  of  wr,  excellent  agreement 

for  these  two  moments  indicates  that  the  basis  of  the  present  formulation  is  sound.  It  is 

worth  noting  that  a  naive,  Eulerian-based  estimation  of  the  average  collision  velocity 

4      ,  , 
based  on  the  incoming  particle  flux,  t  TR  nj,  and  the  available  collision  area,  27tR  , 

2  71 

would  lead  an  erroneous  value  of  ~  TR  instead  of  77:  TR  for  <-w,>. 


Table  3.1        Comparison  of  the  collision  velocity  among  monosized  particles 
between  prediction  and  direct  numerical  simulation  in  laminar 
shear  flow. 


Particle 
radius 

<-wr> 
simulation 

<-wr> 
prediction 

<w2  > 

r 

simulation 

<w2  > 
r 

prediction 

0.009 

5.6421X10-3 

5.6549x10-3 

3.7018X10-5 

3.7029x10-5 

0.0125 

7.8638X10-3 

7.8540x10-3 

7.1529X10-5 

7.1429x10-5 

0.02 

1.2425X10-2 

1.2566x10-3 

1.8012X10-4 

1.8286xl0"4 
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In  a  laminar  shear  flow,  wr=-uxcos9  and  the  collision  angle  9C  between  the 
collision  velocity  w  and  the  axis  connecting  the  center  of  colliding  particle  pair, 
Tt/2<ec<Tt,  happens  to  equal  to  0  as  shown  in  figure  3.1.  In  the  simulation,  0C  can  be 
easily  evaluated  so  that  the  pdf  for  9C  is  readily  obtained.  Since  (-wr)       sin0d0d<i)  is 

wr<0 

proportional  to  the  probability  of  the  particle  striking  the  target  surface  at  (0,  (j>)  with  a 
solid  angle  dco=d0d(j).  it  is  readily  seen  that 

(-wr)dA  =  -rR3sin20cos0  cos(|>  d0d<|>  for  (0,  <|>)eQ, 

=  0  elsewhere  (3.35) 

where  Q  is  the  region  defined  by  (|  <0<tc,  -|  <(j)<| )  and  (O<0<| ,  |  <<)><y  )  in  which 

wr<0.  Interpreting  sin20cos0sin<(>  as  proportional  to  the  probability  density  function  of 
incoming  particles  striking  the  target  particle  near  (0,  <|>),  further  integrating  over  <j>  from  0 
to  27i,  treating  0<9<tc/2  as  the  same  as  n/2<Q<n  due  to  symmetry,  and  using  appropriate 
normalization,  one  obtains  the  pdf  for  the  collision  angle  0c 


p(0c)  =  -B  sin^0ccos0c  for  tc<0<tc/2 

=  0  for  O<0<7t/2  (3.36) 

where  B  is  a  normalizing  factor  which  is  3  for  0C  in  radian  and  ti/60  for  0C  in  degree.  The 

average  angle  is  <0C>  =j  0cpdf(0c)d0c  =  2.237  radian  or  128.17°.     The  numerical 

simulation  gives  <0c>  =  128.13°.  Figure  3.4  compares  the  pdfs  based  on  the  above 
prediction  and  the  direct  numerical  simulation  using  ten  thousand  particles  of  radius 
0.009  over  40  realizations  of  initial  particle  positions.  Excellent  agreement  is  observed 
for  the  whole  region  7i/2<0c<7t.  The  present  theoretical  framework  is  thus  validated  in  the 
laminar  shear  flow. 
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Figure  3.4  Comparison  of  the  pdf  s  for  the  collision  angles  0C  between 
the  theory  and  the  simulation  in  a  laminar  shear  flow. 


3.4.2    Collision  Rate  and  Velocity  in  a  Gaussian,  Isotropic  Turbulence 


For  particle  collision  in  isotropic  turbulence  which  is  represented  using  equation 
(2.6),  the  dimensionless  collision  kernel 


au  '  R3(e/v),/2  (s/v),/2 


(3.37) 


is  evaluated  using  NR=500  realizations.  The  dissipation  rate  s  is  based  on  the  ensemble 
average  of  the  same  NR  realizations.  In  each  realization,  ex  tt{Q  and  e/v  are  computed  in 

the  flow  field  with  654  statistically  independent  data.  It  is  found  that  au  varies  slightly 
with  the  number  of  the  Fourier  modes,  Nk,  used  in  the  turbulence  representation. 
Furthermore,  it  weakly  depends  on  the  power  u  in  equation  (2.10)  with  which  the  random 
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wavenumber  is  chosen.  Figure  3.5  shows  the  variation  of  an  as  Nk  increases  from  100 

to  4096  for  u=1.2  and  1.4  for  r\  0=  0.05  (Rex=53.8).  The  result  seems  to  converge  after 

*  * 

Nk=2048;  an(u=1.2)  converges  to  1.2935  while  an(u=1.4)  converges  to  1.2945  with  a 

difference  less  than  0.1%.  They  are  both  very  close  to  1.2944,  the  value  predicted  by  ST. 
It  is  worth  commenting  that  due  to  the  rapid  decrease  of  the  energy  spectrum  E(k)  in  the 
high  k  range,  not  all  wavenumbers  that  are  randomly  generated  need  to  be  used  in 
equation  (2.10).  For  those  high  wavenumber  modes  that  correspond  to  very  small 
amplitudes,  they  can  be  simply  neglected  without  compromising  the  accuracy.  Thus  the 
effective  number  of  the  random  modes  is  smaller  than  the  specified  N^.  From  the 
simulations,  the  effective  number  of  the  modes  for  u=1.4  is  roughly  three  times  that  for 
u=1.2  for  the  same  N^.  Hence  the  fact  that  it  takes  larger  values  of  Nk  for  u=1.2  than  for 
u=1.4  to  get  a  converged  result  is  not  surprising.  It  is  also  interesting  to  note  that,  using 

Nk=1024  and  u=1.2,  a  higher  value  of  Re^  (=124.4  using  r\  0=0.01)  gives  a*,  =1.2899  as 
* 

opposed  to   a,,  =1.2912  at  Rex=53.8.    This  clearly  shows  that  shape  of  the  energy 

* 

spectrum  E(k)  in  isotropic  turbulence  has  little  effect  on  a,,  .  With  Nk=200  and  u=1.2, 
* 

a,,  =1.2756  which  is  only  1.45%  less  than  the  converged  result.  Thus  in  the  following 
direct  numerical  simulation  for  particle  collision,  Nk=200  and  u=1.2  will  be  used. 

To  validate  the  present  analyses,  direct  numerical  simulations  of  particle  collision 
in  isotropic  turbulence  are  performed.  The  turbulence  in  the  periodic  box  of  length  n  has 
the  following  characteristics:  root-mean-squared  turbulence  velocity  uo=1.0,  integral 
length  scale  Ln=0.594,  Taylor  micro  length  scale  X.=0.277,  Kolmogorov  length  scale 
11=0.023,  Kolmogorov  time  scale  Tk=(v/s),/2=0.073,  and  Rej=40.1.  With  the  modified 
Fourier  modes  in  the  periodic  box,  the  dimensionless  collision  kernel  based  on  the 
prediction  is  a,,  =1.2711  as  opposed  to  a,,  =1.2756  without  the  modification  of  the 
wavenumbers.  Initially,  five  hundred  particles  of  radius  r=0.005  (r/n=0.217)  are 
randomly  distributed  in  the  box  of  volume  n3  with  a  particle  volume  concentration  of 
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8.44x1 0"6.  These  inertialess  particles  are  advanced  using  local  fluid  velocity.  A  total 

time  period  of  T=6  (i.e.  T/xk=82.2)  is  simulated  with  time  step  At=0.01.  The  average 

number  of  collisions  within  the  entire  period  of  the  simulation  0<t<T=6  in  each 

realization  is  about  0.4271  <1  based  on  our  new  expression  of  the  collision  rate  in  §6.2. 

To  obtain  a  reliable  statistics,  a  total  of  10,000  realizations  is  used  to  generate  an 

* 

estimated  collision  pairs  Nc=4271.  Figure  3.6  (a)  shows  the  simulation  result  for  a,,  in 
0<t/xk<82.2  and  the  predicted  value  of  1.2711.  There  is  a  significant  noise  in  the 
numerical  simulation  data  due  to  the  small  number  of  collisions  within  each  time  step 

despite  the  fact  that  10,000  realizations  have  been  used.  A  cumulative  average  is  applied 

* 

to  reduce  the  noise.  Figure  3.6  (b)  shows  the  cumulative  time  average  of  an  from  the 

simulation.  The  direct  numerical  simulation  is  less  than  the  predicted  value  (1.271 1)  by 

about  2.5%.  This  difference  may  be  caused  by  a  small  effect  of  particle  size,  since  the 

theoretical  predictions  are  expected  to  be  accurate  only  in  the  small  size  limit,  i.e.  r/r|«l, 

while  the  particle  size  in  numerical  simulation  is  r/n=0.217.  The  effect  of  finite  particle 

* 

size  on  the  collision  rate  will  be  explored  in  Chapter  6.  The  agreement  for  otn  between 
the  simulation  and  the  prediction  is  quite  good. 

Table  3.2  shows  the  results  of  the  prediction  and  the  numerical  simulation  for  the 
two  moments  of  particle  collision  velocity,  <-w,>  and  <wr  >,  in  isotropic  turbulence.  The 
predicted  values  are  based  on  u=1.2  and  Nk=200  in  the  periodic  box  with  NR=400 
realizations.  The  particle  sizes  are  r/n  =  0.2175  and  0.4875,  respectively.  Satisfactory 
agreements  are  obtained.  On  the  other  hand,  ST  (1956)  gave  an  expression,  equation 
(3.9),  for  the  Eulerian-based  average  particle  radial  collision  velocity.    Using  that 

<-W,> 

expression,  one  would  obtain  =  0.206.  The  present  prediction  for  the  periodic 

domain  with  NR=400,  after  ensemble  average,  gives  a  value  of  0.328  for  this 
dimensionless  radial  collision  velocity  which  agrees  rather  well  with  the  direct  numerical 
simulation  show  in  Table  2.  Since  the  direct  numerical  simulation  is  based  on  the 
tracking  of  particle  trajectories,  the  probability-based  prediction  for  <-w,>  is  consistent 
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with  the  Lagrangian  approach.  The  Eulerian-based  approach  for  <-w,>  does  not  take  into 
account  the  bias  of  the  collision  probability  toward  the  high  values  of  -wr;  hence  it  gives  a 
lower  value  for  <-w,>. 


1.0E+2  1.0E+3  1.0E+4 


Nk 

Figure  3.5       Convergence  ofan  as  ->8. 


Table  3.2  Comparison  of  the  collision  velocity  among  monosized  particles 
between  prediction  and  direct  numerical  simulation  in  Gaussian, 
isotropic  turbulence. 


Particle 
radius  (r/r|) 

<-wr> 
simulation 

<-wr> 
prediction 

<w2  > 

r 

simulation 

2 

<w  > 
r 

prediction 

0.217 

4.5430X10-2 

4.5611x10-2 

2.7356x10-3 

2.7296X10-3 

0.4875 

1. 0093X10-1 

1.0040x10-" 

1.3405x10-2 

1.3224x10-2 
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Figure  3.6  Particle  collision  kernel  an  from  direct  numerical  simulation  in  Isotropic 
turbulence,  (a)  Variation  of  the  ensemble  average,  (b)  Cumulative  time  average. 
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3.4.3    Collision  in  a  Rapidly  Sheared  Homogeneous  Turbulence 

The  statistics  of  the  turbulence,  such  as  four  non-zero  components  of  the 
Reynolds  stress  tensor,  generated  using  RDT  and  the  random  Fourier  modes  was  first 
checked  against  the  analytical  values  at  short  times  (Tt«l)  given  by  RDT  (Townsend, 
1976,  p. 84)  to  ensure  the  correct  implementations.  It  was  found  through  simulation  that 
the  development  of  the  sheared  turbulence  with  time  does  not  depend  on  the  Reynolds 
number,  Re^,  of  the  initial  state  or  the  energy  spectrum  at  t=0.  The  turbulence  part  of  the 
flow  in  the  RDT  is  determined  by  the  total  shear  S=rt.  The  complete  flow  field,  sheared 
turbulence  plus  the  mean  shear  flow,  thus  depends  on  two  parameters:  T  and  S=rt.  Since 
the  collision  rate  depends  on  the  spatial  structure  of  the  turbulence  while  the  turbulence 
structure  becomes  increasingly  anisotropic  as  the  total  shear  S  increases,  finding  an 
appropriate  representation  of  the  collision  rate  in  the  homogeneous  turbulent  shear  flow  is 
not  trivial. 

To  develop  a  better  understanding  on  the  dependence  of  the  collision  rate  on  the 
turbulent  flow,  it  is  instructive  to  examine  the  collision  rate  due  to  the  sheared  turbulence 
but  without  the  contribution  from  the  mean  shear  T;  this  part  of  the  collision  kernel  will 
be  denoted  as  an.    Another  word,  an  is  evaluated  using  the  flow  field  given  by 

equation  (2.18)  but  without  the  last  term  rx28n  which  is  deterministic.  Following  the 
* 

definition  for  an  in  the  isotropic  turbulence  case,  we  use  the  dissipation  rate  (e/v)1/2, 
which  increases  with  S,  as  the  turbulence  characteristic  shear  rate  to  normalize  an, 

d|  i 

a' 1  ~  R3(e/v)l/2     (without  contribution  from  mean  shear  r).  (3.38) 
_  * 

It  is  easily  seen  that  d, ,  depends  on  the  total  shear  S=rt  only.  Figure  3.7  (a)  shows  the 

.  .         .*  * 
variation  of  an   on  S.   In  the  high  S  limit,  the  simulation  results  suggest  that  dn 

approach  a  constant  of  1.1676  in  the  form  of 
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a,,  -1.1676 +  0.659/S        forS>10.  (3.39) 

In  this  rage  of  0<S<50  turbulence  goes  from  isotropic  to  strongly  anisotropic  state.  The 
turbulence  intensity  becomes  dominated  by  the  longitudinal  component  in  the  x\- 
direction.  Based  on  RDT  prediction  which  neglects  the  nonlinear  interaction  among 
different  scales,  the  longitudinal  component  of  the  turbulence  makes  up  92.7%  of  the 
total  kinetic  energy  at  S=rt=20.  The  anisotropy  invariant  lib  reaches  a  value  of  -0.25  at 
S=20.  At  such  large  values  of  total  shear,  S,  the  RDT  results  may  become  a  little 
inaccurate  in  comparison  with  DNS  result.  However,  for  the  present  work,  since  the 
RDT-based  turbulence  is  strongly  anisotropic,  it  serves  our  purpose  to  evaluate  the  effect 
of  such  anisotropy  in  the  turbulence  on  particle  collision  rate.  When  the  collision  kernel 

is  normalized  using  (e/v)1/2,  the  maximum  difference  in  dn  is  only  about  15%  at  the 

* 

extreme  limit  of  S->oo.  If  the  isotropic  result  is  used,  dn  -1.2944,  only  a  maximum 

error  of  10%  will  result  from  the  limiting  value  of  d,,  at  very  large  total  shear.  This 

* 

encouraging  result  will  be  used  to  develop  an  approximation  for  a,,  in  the  presence  of 
strong  mean  shear. 

Now  consideration  is  given  to  the  particle  collision  kernel  with  contribution  from 
both  the  turbulence  part  and  the  mean  flow  shear  part  as  the  flow  field  is  given  by  the 

entire  right-hand-side  of  equation  (2.18).  Figure  3.7  (b)  shows  a*,  as  a  function  of  the 
ratio  of  the  mean  shear  rate  to  the  turbulence  shear  rate,  I7(-  )1/2,  for  a  range  of  shear 

rates.  For  a  given  T,  (e/v)1/2  increases  as  t  or  S  increases  so  that  the  abscissa  mainly 
reflect  the  effect  of  t  or  S.  To  see  the  effect  of  mean  flow  shear  rate,  T  is  made 
dimensionless  by  the  turbulence  mean  shear  rate  at  t=0  in  the  isotropic  state  (8o/v)1/2. 
Four  values  of  r/(eo/v)1/2  are  used:  0.87,  2.18, 4.36,  and  8.72;  they  correspond  to  the  solid 
symbols  in  figure  3.7  (b).  The  open  symbols  specifically  correspond  to  the  cases  when 
S=0  (or  t=0)  so  that  the  turbulence  is  still  isotropic.  The  dimensionless  values  of  a*, 


54 


seem  to  collapse  reasonably  well.  At  large  values  of  I7(~  )1/2,  the  curve  becomes  linear 

as  the  collision  rate  is  increasingly  dominated  by  the  mean  shear  so  that  a,,  scales 
linearly  with  T  as  Smoluchowski's  theory  predicts.  A  small  degree  of  scattering  exists  for 

r/(~  )1/2<2.  This  is  caused  by  the  15%  maximum  variation  of  dn  over  the  whole  range 

of  total  shear  S,  as  already  demonstrated  in  figure  3.7  (a).  There  are  two  possibilities  for 

r/(~  )I/2  «1:  i)  both  f  and  S=ft  are  very  small  even  comparing  with  (so/v)1/2  of  the  initial 


isotropic  turbulence;  ii)  T  is  finite  but  S=Tt  is  large  so  that  (s/v),/2  become  large 

1  /9 

comparing  with  (so/v)    .  In  the  first  case,  the  turbulence  is  closer  to  the  isotropic  state 

.*  * 
so  that  a,,  -1.2944.  In  the  second  case,  the  turbulence  is  highly  anisotropic  and  a,, 

* 

->  1.1 676  for  very  large  values  of  S.  No  attempt  is  made  to  correlate  dn  with  the  total 

shear  S  since  S  is  not  relevant  in  practical  applications  such  as  in  turbulent  pipe  flow. 
_  * 

Using  dn  =1.2944  in  the  isotropic  state  to  represent  the  turbulence  part  of  the  collision 

kernel,  a  simple  interpolation  is  proposed  to  represent  the  effect  of  mean  shear  and 
turbulence  on  the  collision  rate, 

(3.40, 

This  interpolation  is  shown  by  the  solid  line  in  figure  3.7  (b).  A  satisfactory  agreement 
can  be  observed.  It  is  also  worth  commenting  that  when  T  is  normalized  using  the  eddy 
turn  over  time,  3uQ  Izq,  of  the  initial  isotropic  turbulence,  the  shear-rate  parameter 

r*  =  r3u^/80,  (3.41) 
corresponds  to  T*=36.98  and  92.45  for  I7(eo/v)l/2  =0.87  and  2.18.  In  a  low  Reynolds 
number  turbulent  channel  flow,  the  maximum  value  of  T*  is  around  35  at  a  location  of  10 
wall  units  away  from  the  wall  (Kim  et  al.  1987).  Hence,  the  parameters  used  in  figure  3.7 

(b)  are  relevant  to  practical  situations.  Although  a  10%  error  exists  for  a*,  given  by 
(3.40),  the  advantage  of  this  simple  approximation  is  that  the  total  shear  S=Tt  is  absent 
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and  it  depends  only  on  the  spatial  structure  of  the  turbulence  characteristics  by  the  scalar 

1  f) 

quantity  (s/v)  . 

Figure  3.8  shows  the  radial  collision  velocity,  <-w,>,  as  a  function  of  17(~  )1/2,  for 

£0  in 

different  values  of        )    .  Similarly,  the  data  seem  to  collapse  approximately  into  one 

curve  in  the  range  studied.  The  following  interpolation  fits  the  radial  collision  velocity 
well, 

<-Wr>         fn  r7t        r  1.7,1/17 

R(8/v)./2  » (0-3245    +[10(£/v),/2]    }     •  (3-41) 

For  the  root-mean  squared  radial  collision  velocity,  <w2  >,  the  behavior  is  very  similar  to 

that  of  <-w,>.  While  the  data  is  not  shown  here,  the  following  interpolation,  which  fits 
the  data  quite  well,  is  provided, 

<  w2>1/2 

 [  /n  ^o,  1.7  ,  r0I338_L_r  ,1.7,1/1.7 

R(s/v)1/2  -  {03681     +[  (e/v)I/2  ]    }      *  (142) 
The  application  of  the  above  expressions  to  turbulent  channel  flow,  which  is 
different  from  the  homogeneous,  rapidly  sheared  turbulence,  is  yet  to  be  tested.  However, 
since  the  rapidly  sheared  turbulence  at  large  total  shear  S  considered  in  this  work  is 
already  drastically  different  from  the  isotropic  turbulence,  there  are  reasons  to  expect  a 

similar  dependence  of  cti  i,  <-w,>,  and  <w2  >  on  T  and  (s/v)1/2  in  a  turbulent  channel  flow 
which  is  of  significant  industrial  importance. 


3.5     Summary  and  Conclusions 

A  theoretical  framework  has  been  developed  to  evaluate  the  collision  rate  and 
collision  velocity  of  small  particles  in  general  turbulent  flows.  The  classical  results  for 
the  collision  kernels  by  Smoluchowski  (1917)  for  laminar  shear  flow  and  by  Saffman  & 
Turner  (1956)  for  isotropic  turbulence  are  recovered.    The  present  approach  differs 
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significantly  from  that  of  Saffrnan  &  Turner  in  that  the  ensemble  average  is  taken  after 
the  collision  rate  for  a  given  flow  realization  is  calculated.  This  avoids  the  assumption  of 
isotropy,  as  needed  in  Saffrnan  &  Turner,  and  allows  for  the  evaluation  of  the  collision 
rate  in  general  turbulent  flows.  The  collision  velocity  predicted  in  both  laminar  shear 
flow  and  isotropic  turbulence  agrees  well  with  the  direct  numerical  simulation  for  particle 
collisions.  The  present  theory  is  used  to  evaluate  the  collision  rate  and  velocity  of  small 
particles  in  the  rapidly  sheared  homogeneous  turbulence.  It  is  found  that  the  effect  of  the 
turbulence  structure  on  the  collision  rate  and  velocity  can  be  mainly  captured  by  using  the 
mean  turbulent  shear  rate  (e/v)  .  The  combined  effects  of  the  mean  shear  and  the 
sheared  turbulence  on  the  collisions  are  elucidated.  Simple  interpolations  are  developed 
to  evaluate  the  collision  rate  and  velocity  for  arbitrary  mean  flow  shear  rate  and 
turbulence  mean  shear  rate. 
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Figure  3.7  (a)  Normalized  particle  collision  kernel  in  a  highly  sheared  homogeneous 
turbulence  without  contribution  from  the  mean  shear  rate. 
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Figure  7  (b)  Collision  kernel  in  a  rapidly  sheared  homogeneous  turbulence. 
eq  denotes  the  dissipation  rate  of  the  isotropic  state  at  t=0. 
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Figure  3.8  Collision  velocity  in  a  rapidly  sheared  homogeneous  turbulence. 
s0  denotes  the  dissipation  rate  of  the  isotropic  state  at  t=0. 


CHAPTER  4 

EFFECT  OF  INERTIA  ON  THE  PARTICLE  COLLISION  VELOCITY  AND 
COLLISION  ANGLE  IN  GAUSSIAN  TURBULENCE 


Numerical  simulation  of  particle  collisions  in  isotropic,  Gaussian  turbulence  is  carried  out 
to  study  the  effects  of  particle  inertia  on  the  turbulence  induced  collision  velocity  and 
collision  angle.  An  asymptotic  analysis  is  developed  to  predict  the  effect  of  the  small 
particle  inertia  on  the  collision  velocity.  A  PDF-based  method  is  also  implemented  to 
predict  the  effect  of  the  small  particle  inertia  on  the  collision  velocity.   The  analysis 

predicts  that  the  increase  in  the  collision  velocity  is  proportional  to  the  terms  (Pxk)~' 

and  (Pxk)  ,  in  which  P  is  the  particle  response  time  and  xk  is  the  Komogorov  time 
scale.  Good  agreement  is  obtained  among  three  different  methods:  the  numerical 
simulation,  asymptotic  prediction,  and  the  PDF-based  method,  for  the  collision  velocity 
of  particles  with  small  inertia.  The  effects  of  both  particle  inertia  and  size  on  the  collision 
angle  have  also  been  investigated.  It  is  found  that  collision  angle  can  be  significantly 
affected  by  both  inertia  and  size.  Average  collision  angle  and  the  variation  of  collision 
angle  increase  as  the  inertia  increases,  while  the  decrease  in  the  size  will  also  reduce  both 
the  average  collision  angle  and  the  variation  of  collision  angle. 


4.1  Introduction 


Particle  collisions  can  take  place  through  a  variety  of  collision  mechanisms. 
Particles  may  be  brought  together  by  inertial  centrifugal  and  gravitational  forces.  The 
velocity  gradient  is  responsible  for  collisions  among  particles  in  laminar  shear  flow  and  is 
a  dominant  factor  for  particle  collisions  in  turbulence.  Brownian  motion  and  Van  der 
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Waal  forces  can  also  affect  the  collisions.  One  of  the  common  characteristics  of  these 
collision  mechanisms  is  the  ability  to  generate  a  negative  relative  velocity  between  two 
particles  along  the  line  connected  with  two  particle  centers.  Of  significance  are  the 
statistics  of  the  relative  velocity  of  two  particles  at  contact  (hereinafter  referred  as 
collision  velocity)  and  the  relative  angle  between  the  velocities  of  two  particles  at  contact 
(hereinafter  referred  as  collision  angle).  For  solid  particles,  collision  causes  momentum 
exchange  among  particles  and  may  influence  the  local  and  global  particulate  flow 
conditions.  The  suspended  solid  particles  are  produced  during  the  manufacturing  process 
of  materials,  the  combustion  of  heavy  fuel  and  pulverized  coal.  The  collisions  among 
suspended  particles  will  greatly  affect  the  properties  of  particulate  flow  and  thus  influence 
the  performance  of  the  systems.  Adhesion  may  occur  when  two  colliding  solid  particles 
have  small  collision  velocity,  while  large  collision  velocity  will  produce  large  restitutive 
rebounce  velocities  for  the  colliding  particles  and  may  overcome  the  bonding  force 
determined  by  cohesive  surface  energy.  The  collision  also  affects  the  angular  momentum 
of  particles.  Tangential  velocity  components  will  be  imposed  to  the  colliding  particles  if 
the  collision  angle  012  is  not  equal  to  zero.  The  particle  rotation  will  be  either  enhanced 
or  reduced  depending  on  the  collision  angle  012,  which  determines  the  orientation  of  the 
particle  tangential  velocities  upon  collision.  Gunaraj  et  al.  (1997)  peroformed 
experiments  measuring  the  collision  angle  6]2  between  heavy  particles  (silica  sand 
particles  with  1000-1500  um  in  diameter)  in  turbulent  channel  flow. 

The  collision  velocity  is  also  a  critical  quantity  in  determining  the  collision  kernel 
and  collision  rate.  In  uniform  shear  flow,  Smolochowski  (1917)  integrated  the  volume 
flow  rate  using  collision  velocity  along  the  spherical  surface  of  collision  radius  R 
(=  Tj  +  rj ,  where  t{  and  rj  are  the  radii  of  collision  pair  respectively)  and  obtained 
collision  kernel  for  uniform  shear  flow.  In  turbulent  flow,  particle  collisions  are  affected 
by  local  turbulence  fluctuations.  In  the  absence  of  any  other  force  creating  relative 
motion  between  particles,  these  local  fluctuations  entirely  control  relative  velocity 
distribution  which  determines  the  collision  kernel.  In  most  of  the  previous  theories,  the 
collision  kernel  in  a  particle-laden  system  is  linearly  related  to  the  collision  velocity,  the 
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first  moment  of  the  relative  velocity  distribution  of  colliding  particles.  Thus  the  scale  and 
form  of  these  particle  relative  velocity  distributions  directly  determine  the  collision  rate. 
The  fact  that  the  collision  rate  and  collision  kernel  in  previous  studies  are  all  proportional 
to  the  collision  velocity  explains  the  critical  importance  of  the  prediction  of  collision 
velocity  (Saffman  &  Turner  1956,  Williams  &  Crane  1983,  Kruis  &  Kusters  1996,  et  al). 

The  particle  collision  has  been  studied  in  two  different  ways.  One  is  based  on 
Eulerian  method,  which  is  favorable  to  theoretical  prediction  and  was  followed  by 
previous  studies  on  the  theoretical  evaluation  of  both  collision  velocity  and  collision 
kernel.  The  other  is  to  use  Lagrangian  tracking  in  direct  numerical  simulation  of  particle 
collisions.  For  inertialess  particles  in  both  uniform  laminar  shear  flow  and  isotropic 
turbulence,  the  collection  of  collision  velocities  using  Lagrangian  method  shows  a  biased 
average  in  favor  of  those  regions  where  wr  is  high,  as  shown  clearly  in  Chapter  3. 
Saffman  &  Turner  (1956)  obtained  the  collision  velocity  using  average  integration  on  the 
spherical  collision  surface  by  Eulerian  method  without  considering  the  influence  of  the 
biased  average.  Although  their  formulation  correctly  predicts  the  velocity  correlation 
between  two  points  for  turbulent  flow,  it  has  been  incorrectly  interpreted  as  the  collision 
velocity  upon  collisions.  The  effect  of  particle  inertia  on  the  collision  velocity  was 
recently  studied  also  by  Sundaram  and  Collins  (1997).  They  used  two  statistical 
variables:  g(R),  the  particle  radial  distribution  function  at  contact,  and  P(w/R),  the 
conditional  relative  velocity  probability  density  function  at  contact,  where  w  is  the 
relative  velocity  of  two  particles  at  contact.  It  is  noted  that  in  Sundaram  and  Collins' 
work,  in  additional  to  the  assumption  of  the  exponential  form  for  P(w/R),  the  g(R)  was 
evaluated  based  on  their  numerical  extrapolation  form  for  distances  smaller  than  the  first 
numerically  evaluated  point.  However,  the  usefulness  of  g(R)  and  P(w/R)  is  rather 
questionable. 

In  this  work,  we  study  the  effect  of  particle  inertia  on  increasing  particle  collision 
velocity  and  collision  angle  by  carrying  out  numerical  simulation  of  particle  collisions  in 
an  isotropic,  Gaussian  turbulence,  and  developing  an  asymptotic  analysis  to  predict  the 
effect  of  small  particle  inertia  on  the  collision  velocity  based  on  Lagrangian  description. 
An  integration  method  is  also  implemented  to  predict  the  effect  of  the  small  particle 
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inertia  on  the  collision  velocity.  The  increase  in  the  particle  collision  velocity  due  to 
inertia  predicted  by  both  the  asymptotic  analysis  and  integration  method  is  in  agreement 
with  the  numerical  simulation  results.   The  analysis  predicts  that  the  increase  in  the 

collision  velocity  contains  terms  that  are  proportional  to  (Pxk)_1  and  (pxk)-2.  Good 
agreement  is  obtained  by  the  numerical  simulation,  asymptotic  prediction  and  the 
integration  method  for  the  collision  velocity  of  particles  with  small  inertia.  The  effects  of 
particle  inertia  and  size  on  the  collision  angle  have  also  been  investigated.  It  is  found  that 
collision  angle  are  significantly  affected  by  both  inertia  and  size.  Average  collision  angle 
and  the  variation  of  collision  angle  increase  as  the  inertia  increases,  while  the  decease  in 
the  size  will  reduce  both  the  average  collision  angle  and  the  variation  of  collision  angle. 


4.2  Saffman  &  Turner's  Analysis 

In  Saffman  &  Turner's  (1956)  prediction  of  the  collision  rate  among  inertialess 
particles,  the  volume  flow  rate  of  particles  towards  a  target  particle  was  evaluated  by  an 

integration,  <JWr<0  -wrdS>,  where  wr  is  the  radial  component  of  the  relative  velocity. 

The  integration  was  carried  out  over  the  entire  spherical  surface  of  radius  R,  as  shown  in 
figure  4. 1 .  Because  of  the  continuity  of  fluid  flow, 

<Jwr<0  -WrdS>  =  <lentire  sphere  K|dS> 

=  lentire  sphere  <KI>dS.  (4.1) 

For  isotropic  turbulence,  which  is  true  at  small  scales  for  anisotropic  turbulence,  the 
above  becomes 

<Jwr<o  -wrdS>=  27tR2  <|wx|>  (4.2) 

where  x  is  the  local  coordinate  parallel  to  the  line  connecting  the  centers  of  the  colliding 
particles  as  shown  in  figure  4.2  and  <  >  denotes  ensemble  averaging. 
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a  y 


Figure  4. 1   Sketch  of  particle  collision  scheme 


Figure  4.2  Sketch  of  particle  collision  velocity  and  collision  angles  (6 12  &  0C) 
The  separation  distance  is  x2  -  x(  =  ±R  =  ±(rj  +  r2) . 
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In  the  limit  of  zero  particle  inertia,  particle  velocity  v  follows  the  fluid  velocity  u 
very  well  so  that  the  radial  component  of  relative  velocity  w  =    -  v2  upon  collision  is 

|wxHulx-u2x|~R]^|=  RCt^--)172  (4-3) 
ox         15tx  v 

and  <w?>  =  5  R2  <  (— )2  >  =  —  R2-  (4.4) 

dx  15  v 

where  wx  is  the  collision  velocity. 

For  particles  with  inertia,  the  volume  flow  rate  of  particles  to  a  target  particle  was 
expressed  as  JJJentire  sphere  I  wl  P(w)dw,  where  P(w)  is  the  probability  density  function 

(PDF)  of  w  and  contains  the  effects  of  the  spatial  variations  of  velocity  in  the  fluid,  the 
turbulent  accelerations  due  to  the  relative  motion  with  the  turbulence,  and  gravity.  The 
volume  flow  rate  can  be  integrated  assuming  P(w)  to  be  a  normal  distribution  function. 
For  monodisperse  particles, 

Entire  spherel-|P(-)^-|R3(^),/2.  (4.5) 

Compared  with  the  equation  (4.2),  the  equivalent  collision  velocity  based  on  equation 
(4.5)  can  be  evaluated 

Kl  ~  ^R3  (^),/2/(2nR2)=  ^R(-^),/2  (4.6) 
3         v  3  2nv 

and  <w2>~  -R2-.  (4.7) 

9  v 

It  is  noted  that  equation  (4.7)  is  independent  of  the  inertia  parameter  (5.  For  the  small 
inertia  limit  [5-*oo,  comparing  equations  (4.6)  &  (4.7)  with  equations  (4.3)  &  (4.4),  an 
inconsistency  that  (4.6)  &  (4.7)  can  not  recover  (4.3)  &  (4.4)  occurs,  which  is  physically 
incorrect. 
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4.3  Present  Analysis 

4.3.1  Asymptotic  Analysis 

To  understand  and  predict  the  effect  of  inertia  on  the  collision  velocity  wx ,  it  is 
assumed  that  the  particle  size  is  much  less  than  the  Taylor  micro-length  scale,  and  the 
gravitational  force  is  negligible.  For  simplicity  it  is  assumed  that  turbulence  is  isotropic 
and  Gaussian.  The  particle  dynamic  equation  in  the  absence  of  (or  weak )  gravity  is 

^  -  m  -  v).  (4.8) 

The  above  can  be  written  as 

v  =  u  .  (4.9) 

p  dt  v  ' 

For  large  values  of  p,  repeatedly  substituting  the  left  hand  side  into  the  right  hand  side 
gives 

1  du     1  d2v_       lDn    1,      „  „   ,    1  d2v 


V  =  U  +  — r-  -=  U  +  — (U-V)  VU+ 


p  dt  '  p2  dt2        P  Dt    p  p2  dt2 

~«.I£+0(-L).  (4.10) 
P  Dt        p2  v  ' 

The  variance  of  wx  =  vlx  -  v2x  in  which  particle  1  and  particle  2  are  separated  by  a 

distance  of  |  x2  -  X]  |  =  R  is 

,    ,  1  Dulx     1  Du2x 

var(wv )  ~  var(ulx  -  u->„  —  +  —  ) 

x/  lx      2x  p,    Dt      p2  Dt 


2  |  ,  1  Dulx     1  Du2x  2 

-<(ulx-u2x)  >|X2_X)=±R  +<(--—-—__)  >|X2_X1=±R 


7<(u     -U     V  1  DU'x       1  DU2xw, 

2<(U1X     U2xXpi    ^  -^-5r)>lx2-x1=±R 

=  term  I  +  term  II  +  term  III.  (4.11) 
The  first  term  can  be  expanded  as 
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terail  =  <(ulx-u2x)  >lx2-x,=±R  =<u?x>  +  <u2x>-2<uIxu2x>|X2_Xl=±R 

=  2<  u2  >  -  2<  u2  >f(r=R)  =  2  u2  [l-f(R)]  (4. 12) 

where  u2,  is  the  mean  square  intensity  of  turbulence  and  f(r)  is  the  standard  longitudinal 
velocity  correlation  coefficient.  If  the  particle  size  is  much  less  than  the  Taylor  micro- 
length  scale  X,  the  above  can  be  simplified  to 

term  I  =  <(ulx  - u2x)2  >|  X2  _X|  =±R  ~  2  u2  -  A2  =  -U2  £ 

1     1  2  1         15  v 

for  R«A..  (4.13) 

Term  II  in  equation  (4.1 1)  can  be  similarly  handled, 


kp,    Dt      P2    Dt        lx2-x,-±R      p2    I  Dt  )         p2     I  Dt  ) 


_L<Du,x  Du2x 

p,p2       Dt     Dt     1X2   X'_±R  K  } 


where  <  Du'*  D"2*  >|  R  is  the  spatial  correlation  of        in  the  x-direction  with  a 

Dt    Dt       2    1  K  Dt 

separation  distance  of  R.  Noting  that  var(^5-)  =  var(  P^2x  )  and  expanding  the 
correlation  for  small  R,  we  obtain 

<^>|x    x.±r.<(^>[1+^+.  ]  (4  15) 

where  fD  is  the  correlation  coefficient  of        .  Defining  the  micro  length  scale  of  the 

Dt 

fluid  acceleration  as 
Dt 


1  ^  152fD(0) 
A.2D      2  5x2 

term  II  can  be  expressed  as 


(4.16) 


Dtj     Vp2      p2^    P|p2    (Dt)     U  ^ 
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^O^^T^^F^^       for  R<<^-  (4J7) 

Dt         P,     P2        P,P2       01  kD 
Finally,  equations  (4.13)  &  (4.17)  lead  to 

1    2  e      ,DuXx2     1      1  2      2     .Du^  R2 

IS  v+  w    pT  "  FT  'ST**'  ^ 

-^«— *ir£-ir£>-  (418> 

For  the  like  particles,  Pj  =  p2  =  P , 


lp2£     2  vDulx    Du2  2     Du   2  R2 

—  R  <lu,.-u,.)(  )>+^r<(  )  > -^r~ 

15     v     p  2x     Dt       Dt         p2    V  Dt  X2D 


(i)  (ii)  (iii) 

(4.19) 

where  term  (i)  is  the  mean-squared  collision  velocity  for  inertialess  particles,  terms  (ii)  & 
(iii)  are  contributed  by  the  inertia.  The  net  increase  of  the  mean-squared  collision 
velocity  due  to  the  inertia  can  be  expressed  as 

,2 

Dt       Dt  '      p2  "v  Dt  '  '  X2D 
(ii)  (iii) 


2 

i  2    ,  wDu,„    Du2x  .       2    ,Du„.2    R  . . 

<AWX>  ~    --  <(UK-U2J(^-— ^)>+-y<(— ^) 2>-y-.  (4.20) 


It  is  noted  that  term  (ii)  represents  the  first-order  effect  of  the  inertia  and  term  (iii) 
contains  the  second-order  effect  of  the  inertia.  The  present  prediction,  equations  (4.19)  & 
(4.20),  differs  from  the  equation  (4.7)  in  two  aspects.  First,  the  collision  velocity  in  the 
limit  p~>oo  is  consistent  with  the  equation  (4.4).  Second,  the  present  analysis  includes 
two  additional  terms  in  the  effect  of  inertia,  as  seen  from  equations  (4.19)  &  (4.20)  when 
P]  =  P2 .  The  effect  of  inertia  does  not  vanish  when  P|  =  P2  while  it  does  in  Saffman  & 

Turner's  analysis.  It  should  be  mentioned  that  <Awx>  in  (4.20)  is  computed  based  on 
Lagrangian  tracking  of  all  colliding  particles. 
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4.3.2  PDF-Based  Method 


For  inertialess  particles,  the  collision  velocity  can  be  evaluated  by  the  direct 
integration  (Mei  &  Hu  1998)  in  both  uniform  shear  flow  and  turbulent  flows. 


JWr<owr(-wr)dS 
<wr>=  ~i  —  (4.21) 


and 


2  Jwr<0Wr(-Wr)dS 

<w?>="7  7  ^~  (4-22) 

Jwr<0(-Wr)dS 

where  (-wr)dS  in  the  numerator  can  be  interpreted  as  proportional  to  the  probability  of  a 
particle  striking  the  target  surface  over  an  area  dS.  The  denominator  is  the  normalizing 
factor  for  the  probability.  The  product,  wr(-wr),  clearly  indicates  that  <w,>  is  a  biased 

average  in  favor  of  these  regions  where  (-wr)  is  high. 

Equations  (4.21)  &  (4.22)  may  be  extended  for  particles  with  inertia.  For  small 
inertia,  as  shown  in  §4.3.1,  the  particle  velocity  to  the  first  order  p  is 

1  Du 

v  =  „---.  (4.23) 
The  collision  velocity  may  be  expressed  as 

1  Dulx      1  Du2x 

w  =  wx  =  vlx  -  v2  =  ulx  -u2x  —  +  —  (4.24) 

r        x        lx       2x        ix        2x    P)     ^       p2     Dt  *  > 

r  jo     r  1    Du,„       1    Du2v  . 

fWr<o  -wrdS  =  IWr<0  -(ulx-u2x-— -^  +— -^)dS.  (4.25) 
For  monodisperse  particles,  P,  =  p2  =  P ,  equation  (4.25)  may  be  rewritten  as 

lWr<o  -wrdS  =  jWr<0  tu2x-Uix+^(^L-^L)]dS.  (4.26) 
From  the  equations  (4.21),  (4.24)  and  (4.26),  the  collision  velocity  is  then  found 
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<J.,^u2x-ulx+I(^-^dS> 

<«-,>=   f-rP1  =^  .  (4.27) 

<L,<„[u2,-u,l  +  i(^-I^)]dS> 

Similarly,  the  mean-squared  collision  velocity  is 

<w?>  =     k<°  nD'      nD'  .  (4.28) 

The  integrations  in  right  hand  side  of  (4.28)  are  evaluated  using  Eulerian  method  . 


4.4  Results  and  Discussions 


The  direct  numerical  simulations  of  the  particle  collisions  in  Gaussian  turbulence 
is  performed  to  assess  present  asymptotic  analysis  and  to  predict  collision  angle.  Initially, 

a  given  number  of  particles  are  randomly  distributed  in  the  box  of  volume  rc3.  The 
inertial  parameter  Pxk  varies  from  0.022  to  qo.  The  particle  sizes  D/r|=0.435,  0.957, 
1 .739,  2.522,  4.348  are  considered,  and  the  corresponding  initial  number  of  particles  are 
24000,  10000,  5000,  3000  and  1000  respectively.  The  particle  collision  detection  is 
turned  on  until  all  particles  have  reached  dynamic  equilibrium.  The  overlapping  particles 
are  eliminated  to  avoid  the  false  counting  of  the  collisions  before  the  collision  detection 
turned  on.  The  particles  are  advanced  by  time  step  At=0.01.  The  total  time  steps,  varying 
from  200  to  600  after  the  equilibrium  state,  depend  on  the  inertia.  A  total  of  40 
realizations  is  used  for  smooth  statistics. 

Figure  4.3  shows  the  comparisons  of  direct  numerical  simulations,  asymptotic 
analysis  and  PDF-based  method  for  D/r|=0.957.  Very  good  agreement  is  observed  among 
the  three  methods  for  small  inertia  Pxk>3.6,  so  equation  (4.20)  is  validated  for  prediction 
of  small  inertia  effect  on  the  collision  velocity.  However,  large  deviation  is  observed  for 
Pxk<3.6  since  the  asymptotic  analysis  cannot  be  applied  for  particles  with  large  inertia. 
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The  reason  is  that  the  first  order  approximation  (4.23)  is  no  longer  valid  when  the  P 
becomes  small  for  particles  with  large  inertia.  The  PDF  results  of  collision  velocity  can 
be  found  in  figures  4.4  (a)  &  (b).  The  collision  velocities  are  all  negative,  because 
collision  occurs  only  if  the  pair  of  particles  are  approaching  each  other.  For  not  very 
large  inertia,  say  (3xk>l,  due  to  the  strong  correlation  between  particle  velocity  and  fluid 
velocity,  the  PDF  curves  are  sharp  and  narrowed  with  a  small  variance  and  peak  at  almost 
the  same  location,  while  zero  inertia  curve  is  at  highest  level.  The  PDFs  are  sharper  for 
small  particle  size  (D/rj=0.957)  at  Pxk>1.023,  shown  in  figure  4.4  (a),  compared  with 
figure  4.4  (b).  However,  the  PDF  curves  become  flatter  and  the  location  of  the  peak 
shifts  to  larger  |wr|  for  Pxk<0.073  due  to  little  correlation  between  particle  velocity  and 
fluid  velocity  at  large  inertia. 


Figure  4.3  Comparison  for  the  increased  collision  velocity  due  to  the  particle  inertia 
between  numerical  simulation  and  theoretical  prediction,  D/r|=0.957. 
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Figure  4.4  Effect  of  particle  inertia  on  the  probability  density  function  of  the 
particle  collision  velocity  wr .  (a)  D/r|=0.957.  (b)  D/r|=2.522. 
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The  effect  of  particle  inertia  on  collision  angle  among  monodisperse  particles  can 
be  found  in  figures  4.5  (a-c),  4.6  &  4.7.  It  is  seen  that  the  collision  angles  812  &  0C 

strongly  depend  on  the  inertia.  For  not  very  large  inertia,  (kk>1.023,  the  mean  collision 
angle  §12  is  less  than  20°  at  D/r)=0.957,  as  shown  in  figure  4.5  (a)  &  4.6,  and  it  is  less 
than  30°  at  D/r(=2.522,  as  shown  in  figures  4.5  (b)  &  4.6.  However,  for  very  large 
inertia,  fkk<0.073,  the  mean  collision  angle  0i2  is  greater  than  80  0 ,  which  reflects  little 
biased  orientation  of  collisions  due  to  very  weak  correlation  between  the  particles  and 
turbulent  flow.     It  is  noted  that  the  PDF  curves  for  Pxk>  1.023  are  much  sharper  and 

narrower  than  those  for  Pxk<0.073  and  the  variation  of  collision  angle  increases  as  the 
particle  inertia  increases,  as  shown  in  figure  4.7.  Figure  4.5  (c)  shows  the  effect  of  inertia 
on  0C.  It  is  seen  that  the  increase  in  inertia  results  in  the  collision  angle  0C  approaching 

to  symmetrical  distribution  with  the  mean  45  0 . 

The  effect  of  particle  size  on  the  collision  angle  012  can  be  seen  in  figure  4.8, 
which  clearly  shows  that  the  particle  size  has  a  significant  effect  on  the  collision  angle. 
The  average  collision  angle  912  increases  as  the  size  increases.   For  D/r)<1.739,  the 

maximum  012  is  located  at  somewhere  between  0°and  10°,  and  the  average  012  is 
around  10°.  However,  for  larger  particles  with  large  inertia,  shown  in  figure  4.8,  the 
collision  angle  0 12  is  around  20 0 ,  which  is  within  the  range  of  14 0  ~  32 0  obtained  from 
the  experiments  for  heavy  particles  in  a  turbulent  channel  flow  (Gunaraj  et  al.  1997). 
Since  little  information  is  reported  on  the  turbulent  channel  flow,  the  comparison  is  only 
qualitative. 
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Figure  4.5  Effect  of  particle  inertia  on  the  probability  density  function  of  the 
particle  collision  angle  0I2 .  (a)  D/r|=0.957.  (b)  D/r|=2.522. 
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Figure  4.5  (c)  Effect  of  particle  inertia  on  the  probability  density  function  of  the  particle 
collision  angle  0C,  D/r|=0.957. 


Figure  4.6  Effect  of  particle  inertia  on  the  mean  collision  angle  9i2  in  isotropic, 
Gaussian  turbulence. 
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Figure  4.8  Effects  of  particle  inertia  and  size  on  the  probability  density  function  of 
the  particle  collision  angle  912,  Pxk=1.023. 
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4.5  Summary  and  Conclusions 

Numerical  and  analytical  studies  on  particle  collision  velocity  and  angle  in 
isotropic,  Gaussian  turbulence  have  been  carried  out.  The  effects  of  particle  inertia  and 
particle  size  on  the  collision  velocity  and  collision  angle  are  examined. 

An  asymptotic  analysis  has  been  developed  to  predict  the  effect  of  small  particle 
inertia  on  the  particle  collision  velocity.  The  asymptotic  prediction  for  the  increased 
collision  velocity  agrees  well  with  the  results  of  numerical  simulation  and  PDF-based 

method.  The  collision  velocity  is  found  to  contain  terms  that  are  proportional  to  (Pxk  )_1 
and  (Pxk)-2. 

Particle  collision  angle  strongly  depends  on  particle  inertia  and  size.  The  average 
collision  angle  and  the  variation  of  the  collision  angle  increase  as  the  inertia  increases. 
Both  the  average  collision  angle  and  the  variation  of  the  collision  angle  decrease  as  the 
particle  size  decreases. 


CHAPTER  5 

QUANTIFICATION  OF  PARTICLE  CONCENTRATION  NON-UNIFORMITY 

IN  TURBULENCE 


For  heavy  particles  with  finite  inertia,  small  scale  structure  of  turbulence  induces  particle 
preferential  concentration.  Particles  tend  to  collect  in  the  regions  of  high  strain  rate  and 
low  vorticity,  as  has  been  observed  in  many  Lagrangian  simulations.  Quantification  of 
the  extent  of  the  non-uniform  particle  concentration  distribution  has  been  difficult 
because  the  particle  concentration  based  on  the  Lagrangian  statistics  is  highly  noisy.  A 

method  is  developed  to  calculate  reliably  the  variant  of  the  particle  concentration,  <  n'  >, 
in  which  n  is  the  instantaneous  particle  number  concentration  and  the  prime  is  the 

fluctuation  in  a  small  cube  of  volume  (Ax)3,  from  the  Lagrangian  simulations  of  the 
particle  motions  in  an  isotropic  turbulence.  The  concentration  non-uniformity  is 
measured  using  a  new  parameter  Ap  based  on  the  first  two  lower  order  moments  of  n.  It 

is  observed  that  the  non-uniformity  varies  with  Ax,  suggesting  a  fractal  nature  of  the 
particle  concentration  distribution.  For  particles  with  zero  size,  the  non-uniformity  scales 

with  Ax  as  Ap~  (Ax)_b  (b  around  1)  for  small  values  of  Ax.  However,  for  finite-size 

particles,  the  non-uniformity  reaches  a  maximum  value  (around  order  one  or  less)  at  Ax/D 
=  2.0-3.5,  in  which  D  is  the  particle  diameter;  it  drops  to  zero  as  Ax/D  gets  close  to  one. 

2  1 

It  is  noted  that  the  maximum  value  of  the  <n'  >/<n>  is  attained  when  the  ratio  of 
particle  response  time  to  the  Kolmogorov  time  scale  is  close  to  one,  which  has  been 
previously  observed  for  the  zero-size  case  based  on  other  physical  quantities. 
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5.1  Introduction 

In  many  statistical  problems,  the  lower  order  moments  are  often  of  interest.  For 

large  number  of  monodisperse  particles,  the  collision  rate  is  nc=ann  .    The  first 

moment  of  particle  concentration  is  clearly  the  average  concentration  n ,  i.e.,  the  number 
of  particles  per  unit  volume.  The  second  order  moment  <n  >  can  then  be  expressed  as 

2         -2         ,2       -2        <  n'2  >  ,  ,_ 

<n2>=  n  +<n'z>=n  (1+    _2    )  (5.1) 

n 

where  <  >  denotes  ensemble  average.  The  first  moment  n  does  not  tell  how  the  number 
of  particles  varies  with  the  volumes  at  different  locations  and  gives  no  reflection  of 
turbulence  structure  at  small  scale  levels.  However,  the  second  order  moment,  <n  >, 

consists  of  two  terms:  n  and  <n'  >.  The  variance  of  particle  concentration,  <n'  >, 
measures  the  variation  of  particle  number  density  in  the  cells,  i.e.  the  non-uniformity  of 
the  particle  distribution.  Thus  it  is  fundamentally  important  to  the  particle  collision  rate. 

-  ...  .9 

Since  n  is  a  simple  average  of  particle  concentration  in  a  turbulent  flow  while  <  n'  > 
depends  on  the  turbulence  structure  and  particle  inertia,  etc.,  our  attention  is  given  to  the 

behavior  of  <  n'  >. 

2  _2 

In  this  study,  we  use  the  normalized  variance  of  particle  concentration,  <  n'  >/  n 
to  quantify  the  concentration  non-uniformity  in  turbulent  flows.  The  normalized  variance 
is  obtained  via  direct  numerical  simulation.  As  will  demonstrated  in  §5.2,  the  inherent 

2  -2 

statistical  noise  is  eliminated  from  the  <n'  >/n  .  By  this  new  method  no  assumption 
regarding  the  particle-particle  interactions  is  needed.  The  effects  of  particle  size  and 
inertia  on  the  concentration  non-uniformity  can  be  quantified  separately.  It  is  found  that 
the  concentration  non-uniformity  can  be  greatly  overestimated  if  the  particle  size  is 
ignored.     For  particles  with  zero  size,  the  non-uniformity  scales  with  Ax  as 

<n'  >/<n>  ~  (Ax)    (b  around  1)  for  small  values  of  Ax  for  particles  whose  response 
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time  is  comparable  with  turbulence  Kolmogorov  time  scale,  where  Ax  is  the  length  of  a 
small  cube  of  volume  (Ax)3  used  in  the  Lagrangian  simulations  of  the  particle  motions  in 
turbulence.  However,  for  finite-size  particles,  the  non-uniformity  reaches  a  maximum 
value  (around  order  one  or  less)  at  Ax/D  =  2.0-3.5,  in  which  D  is  the  particle  diameter 
and  r|  is  Kolmogorov  length  scale;  it  drops  to  zero  as  Ax/D  gets  close  to  one.  The 

maximum  value  of  the  <n'2  >/<n>2  is  obtained  when  the  ratio  of  particle  response  time 
to  the  Kolmogorov  time  scale  is  around  one. 


5.2. 1  Inherent  Noise  of  the  Concentration  Non-uniformity 

Suppose  the  computational  domain  L3  is  divided  into  Nc  equal-volume  cells  each 

with  length  Ax  so  that  Nc  =  (L/Ax)3.  Obviously  n  =  N/Nc  is  the  first  moment  of  n. 

For  zero-size,  inertialess  particles,  the  probability  of  a  particle  residing  in  a 
particular  cell  is  p  =  1/NC .  With  N  particles  in  the  domain,  the  probability  of  k  particles 

residing  in  a  specific  cell  of  volume  (Ax)3  obeys  the  Bernoulli  distribution, 


5.2  Quantification  Method 


P(k)  = 


N! 


,N-k 


(5.2) 


k!(N-k)! 


The  second  moment  of  n  is,  by  definition, 


<n2>=  Ik2P(k). 


(5.3) 


k=0 


Thus  the  variance  of  the  particle  concentration  for  a  given  cell  is 


<n,2>  =  <n2>-n2=  £k2P(k)-n2. 


(5-4) 


k=0 


For  N/>(1-  p)»l,  P(k)  given  by  (5.2)  can  be  expressed  as 


(5.5) 
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for  k  in  the  neighborhood  of  Np,  which  is  known  as  DeMoivre-Laplace  theorem.  In  the 
limit  of  large  Np(l-  p),  which  is  very  close  to  n  for  small  p,  an  asymptotic  integration 


gives 


<  n2  >  =  J  k2P(k)  dk  ~  n2  +  n(l-  p) .  (5.6) 
o 


Denoting  <  n'^  >  as  <n,2>  for  zero-size,  inertialess  particles  (P-»<»),  we  obtain 

<n,2>=<n2>-  n2  =  n(l-p).  (5.7) 

2  ~2 

Denoting  r  as  the  ratio  of  <  n'  >  to  n  , 

<n,2> 

r=  —  ,  (5.8) 
n 

it  follows  for  zero-size,  inertialess  particles  that 

_<n'2  >_hZ 

-2 

n  n 


For  n  =0(1)  or  smaller,  the  DeMoivre-Laplace  theorem  does  not  hold  anymore.  In  DNS 
the  number  of  particles  used  is  relatively  small  from  the  point  view  of  smooth  statistics, 

as  mentioned  in  Section  I.  In  DNS,  n  may  be  much  less  than  one  as  cell  size  becomes 
small.  With  finite-size  particles,  the  ultimate  number  distribution  of  particles  is  reduced 
to  either  zero  or  one  per  cell  as  the  cell  size  is  smaller  than  0.5773D,  where  D  is  the 

particle  diameter.  For  n  «0(1),  particle  number  concentration  is  more  likely  to  obey 
Poisson  distribution 

Xke~x 

P(k)=— — ,  (5.10) 
k! 

where  \=Np.  The  second  moment  of  n  is  then 

<n2>=  Zk2P(k)=  Zk2  ,  (5.11) 

k=o  k=o  k! 

and  the  variance  of  the  particle  concentration  for  a  given  cell  is 


81 


<n'2>  =  <n2>-n2=  £k2^^-n2. 

k=o  k! 


ForN»l, 


Ik2  Ik2         i[k(k-i)  +  k]^- 

k=o       k!      k=o       k!      k=0  k! 


=  Zk(k-l)^-+Sk^r 
k=o  k!      k=o  k! 


oo      T>k-2  oo      >  k-1 

=^2e-x  s  J:  +Xe-x  s  *  

kto  (k-2)!  kto(k-l)! 

=  X2  e~x  ek+ke~x  ex 

=  X2+X.  (5.13) 
Since  A.=Np=  n ,  the  variance  of  n  can  be  expressed  as 


n.2>=  l^^-.n2  ~X2+X-n2 
kto  k! 


—2     -  -2  — 

=  n  +  n-n  =  n.  (5.14) 


It  follows  for  zero-size,  inertialess  particles  that 


2 

<  n'oo  >  _    1  ,  « 

roos— 3^  =  ■  (5-15) 

n  n 


2       —  ^ 

For  n  »1,  seen  from  both  equations  (5.9)  and  (5.15),  r^  =<  n'£,  >  In  -»0  so  that 

the  variance  of  the  particle  concentration  distribution  is  very  small  when  it  is  normalized 
by  the  square  of  the  first  moment  of  n.  However,  if  N  is  fixed  while  Nc  takes  a  large 

value  in  the  limit,  we  have  n  «  1  and  r^  ~l/n  »1.  In  a  practical  simulation  of  particle 
motions  in  turbulence,  N  is  always  finite  and  fixed.  Since  we  are  interested  in 
quantifying  the  small-scale  variation  of  the  particle  concentration,  it  is  necessary  to 

reduce  successively  the  size  of  the  cell,  (Ax)  ,  on  which  n  and  <  n'  >  are  computed.  It  is 
noted  that  the  relatively  large  value  of  <  n'^  >  for  inertialess  particles  is  entirely  due  to  the 
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fact  that  there  is  not  enough  particles  in  each  cell  as  the  size  is  reduced  and  r^-l/n 
contains  no  contribution  from  the  preferential  concentration.  As  far  as  the  inertial  effect 

is  concerned,  r~l/n  is  purely  the  inherent  statistical  noise  which  can  be  improved  by 

m  2*1 

using  large  values  of  n  in  principle.  As  the  inertia  of  the  particles  increases,  <  n'  >  will 
contain  an  increasing  effect  of  the  inertia  bias.  The  fact  that  raD~\/n»\  for  inertialess 

2      —  2 

particles  when  n  «  1  implies,  however,  that  r=  <  n'  >/  n   is  dominated  completely  by 

2 

the  inherent  statistical  noise  r^-l/n  .  The  contribution  from  the  inertial  effect  to  <n'  > 
or  r  thus  must  be  isolated  from  the  statistical  noise.  Figure  5.1  shows  the  comparison 

between  the  DNS  results  r^for  inertialess  particles  and  the  inherent  noise  1/n  ,  where 

7  —2 

<n  >/n  at  every  instant  was  recorded  and  a  total  of  300  time  steps  and  40  turbulence 
realizations  are  used  in  our  DNS.  The  details  regarding  turbulence  representation,  the 
characteristics  of  the  turbulence  and  particle  collision  detection  scheme  have  been  given 

9  —2 

in  §2.1.3.  It  clearly  shows  that  <n  >/ n  is  indeed  dominated  by  1/n  which  has  to  be 
subtracted. 

5.2.2  Concentration  Non-uniformity  Parameter 

For  particles  with  finite  inertia,  the  variance  of  the  particle  concentration  can  be 
computed  in  direct  numerical  simulations 
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Figure  5.1  The  normalized  variance  of  the  concentration  as  function  of  cell  size  for 
inertialess  particles  in  isotropic  turbulence.   Pxk  = 1.023. 


The  above  variance  of  the  concentration  non-uniformity  is  dominated  by  the  statistical 
noise,  as  shown  in  figure  5.1.   For  zero-size,  inertialess  particles,   the  above  should 

approach  to  <n'^>=  n(l-  p)  for  n  >0(1),  and  <n'^>=n  for  n  «  0(1).  Defining  the 

normalized  difference  between  above  <n'2(P)>  and  <  n'2  >  as 

n*p  n'jS  (P)  n'l 

<— !->=  < —  >  -  <  — > 

-2  -2  -2 

n  n  n 

np(P)  1-p 

=  <-^T->-  for  n>0(l),  (5.17) 


and 


-2 
n 


np  n'l  (P)  n'l 

< — ->=  < — "  >-< — —> 

-2  -2  -2 

n  n  n 
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n'jS  (P) 
<—  > 

-2 
n 


n 


for  n«0(l). 


(5-18) 


1.20 


080 


0.40 


0.00 


-0.40 


-0.80 


-1.20 


-1.60 


D=0.058.  6q(5.18) 
D=0.032,  Eq(5.18) 
O0.022,  Eq  (5.18) 
zero  size,  Bq  (5.18) 
D=0.058,  Bq  (5.17) 
DO.032,  Bq  (5.17) 
D=0.022,  Bq  (5.17) 
zero  size,  Bq  (5.17) 


10 


Ax  /  T| 


100 


Figure  5.2  The  expressions  (5. 1 7)  &  (5. 1 8)  of  the  concentration  non-uniformity  as 
function  of  cell  size  for  inertialess  particles  in  isotropic  turbulence.  0t  k  =oo. 


The  effect  of  inertia  on  the  variance  of  n  can  be  identified.  Both  equations  (5.17) 

and  (5.18)  represent  the  contribution  to  <n'  >  from  the  inertial  effect  for   n  >  0(1)  and 

n  «  0(1),  respectively.  Figure  5.2  shows  the  results  of  equations  (5.17)  and  (5.18)  for 
isotropic  turbulent  flow.  Equation  (5.17)  obviously  gives  very  small  value  for  large  Ax/D 

for  n  »0(1),  while  equation  (5.18)  on  other  hand  gives  small  values  for  small  Ax/D  for 

n  «  0(1).  For  inertialess  particles,  <n'p>/  n   from  both  expressions  has  0(1)  variation 

2     — 2 

over  cell  size.  Equation  (5.17)  gives  <n'p>/n  ->  -1  as  cell  size  Ax  approaches  to  zero 
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while  equation  (5.18)  gives  <n'p>/n  ->1  as  Ax  increases.  This  is  physically  unsound 

because  the  preferential  concentration  is  an  inertial  phenomenon  and  should  not  appear  in 

turbulent  flows  with  inertialess  particles.  Figure  5.3  shows  the  DNS  results  of  <n'p> 

using  equations  (5.17)  &  (5.18)  in  a  uniform  shear  flow  for  inertialess  particles.  The 
uniform  shear  flow  is  represented  by  u=(Ty,  0,  0)  with  the  velocity  gradient  r=1.0. 

Comparing  with  figure  5.2,  very  similar  behavior  for  <n'p>/n    is  obtained  even  for  a 

uniform  shear  flow.  The  results  are  unphysical  because  there  is  no  mechanism  to  drive 
inertialess  particles  to  form  non-uniform  concentration  in  a  uniform  shear  flow.  Thus, 
one  should  not  expect  to  apply  the  equation  (5.17)  or  (5.18)  to  quantify  the  concentration 
non-uniformity.  However,  the  above  results  gave  us  a  valuable  hint.  That  is,  the  results 
obtained  by  both  equations  (5.17)  &  (5.18)  can  be  interpreted  as  an  additional  statistical 
noise  for  both  inertialess  and  inertial  particles.   It  is  hence  suggested  that  the  actual 

contribution  of  the  inertia  to  <  n'  >  for  an  arbitrary  p  may  be  expressed  as 

<An'j  >  =  [<n,2(p)>-n(l-p)]p  -[<n,2(p=oo)>-n(l-p)]p=a3,  (5.19) 

or 

<A  n'l  >  =  [<n,2(P)>  -  n  ]  p  -  [<n,2(P=oo)>  -  n  ]  p=00 .  (5.20) 

Apparently,  zero  concentration  non-uniformity  is  obtained  for  inertialess  particles  in  both 
turbulence  and  uniform  shear  flows  based  on  equations  (5.19)  &  (5.20).  It  is  noted  that 

<A  n'p  >  is  a  function  of  p  and  n  in  equation  (5.20),  while  <A  n'p  >  is  mainly  determined 
by  p  and  n  in  (5.19).  The  influence  of  p  on  <An'p>  in  (5.19)  may  be  negligible  as  far  as 

a  small  p  is  concerned,  i.e.  <A  n'p  >=  f(P,  n  ),  since  p  appears  in  both  first  and  second 
terms  in  right  side  of  equation  (5.19).  A  large  p,  however,  means  the  large  cell  size  (the 
box  size  is  the  limit),  which  usually  gives  approximately  uniform  particle  distribution.  It 

is  also  noted  that  <An'p>-»0  as  n  ->oo  (or  N-*»).  Thus,  it  is  expected  that  the  effect  of 
inertia  on  the  concentration  non-uniformity  can  be  captured  using  equations  (5.19)  & 
(5.20).   To  quantify  the  extent  of  the  concentration  non-uniformity  using  the  second 
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moment  of  n,  a  dimensionless  quantity  Ap  called  concentration  non-uniformity 
parameter,  is  defined  as 


Ap  =  [— 


<n'l>  Q-^        r<n'2(P  =  oo)>  (1-p) 


]p-[- 


]  B=oo  ' 


(5.21) 


or 


A    -r<n'P>     ln        r<n'2(P  =  °Q)>  1 

Ap  _  L      — =  Jp  -  l  =  Jp=3 


n 


n 


(5.22) 
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Figure  5.3  The  expressions  (5.17)  &  (5.18)  of  the  concentration  non-uniformity  as 
function  of  cell  size  for  inertialess  particles  in  uniform  shear  flow.  p/T=oo. 


Figure  5.4  shows  the  concentration  non-uniformity  parameter  for  the  uniform  shear 
flow  using  equations  (5.21)  &  (5.22).  It  is  seen  that  Ap  is  very  small  for  all  cases  in  a 
uniform  shear  flow.  Figure  5.5  shows  Ap  in  turbulent  flows  for  a  large  range  of  sizes  for 
particles  with  zero  inertia.  It  is  also  shown  that,  for  zero-size  particles,  the  concentration 
non-uniformity  at  small  inertia  (Pxk=7.3)  still  increases  with  the  decreasing  cell  size. 
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The  detailed  discussions  about  the  size  effect  on  the  collision  rate  will  be  given  in  §5.3. 
The  above  results  obtained  by  (5.21)  &  (5.22)  are  reasonable  and  hence  strongly  suggest 
that  equations  (5.21)  &  (5.22)  can  be  used  to  quantify  and  predict  the  particle 
concentration  non-uniformity  induced  by  particle  inertia  in  general  turbulence.  It  should 
also  be  noticed  that  both  equations  (5.21)  and  (5.22)  always  yield  the  same  shape  differed 
by  a  constant  (si).  To  facilitate  the  subsequent  discussions,  equation  (5.22),  which  does 
not  contain  p  and  is  thus  in  a  simpler  form,  is  used  hereinafter. 


Figure  5.4  The  concentration  non-uniformity  parameter  Ap  as  function  of  cell  size  in 
uniform  shear  flow.  r=1.0,  p7T=1.0. 


Figure  5.6  shows  the  effect  of  total  particle  number  N  on  the  concentration  non- 
uniformity  parameter  Ap.  It  is  seen  that  Ap  remains  almost  on  the  same  curve  when  N 

increases  from  25,000  to  50,000  and  100,000  for  all  values  of  Nc(from  23  to  2783). 
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Even  when  N  decrease  to  12,000,  5,000  and  2,000,  Ap  does  not  change  much  from  those 
at  N=  100,000.  This  suggests  that  Ap  is  only  very  weakly  dependent  on  N  and  the  results 
for  Ap  obtained  at  finite  N  (2,000-12,000)  can  be  used  to  represent  Ap  of  N-*».  Hence, 
Ap  defined  by  (5.21)  &  (5.22)  contains  only  the  effect  of  inertia  and  can  be  used 
conveniently  as  a  measure  for  small  scale  concentration  non-uniformity.  It  should  be  also 
pointed  out  that  while  Ap  has  a  minimal  statistical  error  due  to  finite  N,  it  still  contains 
statistical  errors  due  to  finite  number  of  turbulence  realizations. 
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Figure  5.5  The  concentration  non-uniformity  as  function  of  cell  size  for  small  inertia 
particles  in  isotropic  turbulence. 
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Figure  5.6  The  concentration  non-uniformity  as  function  of  cell  size  for  various 
particle  number  densities  at  a  given  inertia  in  isotropic  turbulence.  Px k =1 .023. 


Since  this  chapter  is  focused  on  quantifying  the  concentration  non-uniformity  in 
order  to  investigate  the  effects  of  particle  inertia  and  finite  size  on  the  collision  rate,  the 
"THROW  AWAY"  scheme  is  employed  for  finite-size  particles.  Both  the  tracking  of  the 
motions  of  particles  and  the  collision  detection  are  carried  out  in  a  Lagrangian  frame. 
The  semi-implicit  finite  difference  scheme  with  2nd  order  accuracy  is  implemented  to 
integrate  the  equation  of  particle  motion  and  to  detect  collision. 

To  examine  the  effects  of  the  inertia  and  size  on  particle  concentration  non- 
uniformity  induced  by  the  small-scale  structure  of  turbulence,  N  particles  are  randomly 

seeded  in  a  box  of  n 3  in  an  isotropic  turbulence.  Initially,  the  velocities  of  the  particles 
are  set  to  be  the  same  as  that  of  the  turbulence  at  corresponding  locations.  There  are  two 
categories  of  particles  considered:  zero-size  and  finite-size  particles.    For  zero-size 


5.3    Results  and  Discussions 
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particles,  following  previous  studies  on  the  preferential  concentration,  no  particle 
collision  detection  is  needed  so  that  particles  can  overlap  in  space,  which  is  equivalent  to 
the  use  of  the  "OVERLAPPING"  scheme.  For  finite-size  particles,  however,  collision 
between  particles  has  to  be  detected,  and  the  collided  pair  will  disappear  from  their 
current  size  groups  as  far  as  geometrical  collision  is  concerned.  Periodical  boundary 
conditions  are  imposed  on  both  the  particles  and  turbulence.  Instantaneous  particle 
concentration  fields  are  obtained  by  tracking  2000,  5000,  12000,  25,000  to  50,000  and 
100,000  particles  respectively.  The  box  is  divided  into  Nc  equal-volume  cells  each  with 

length  Ax  so  that  Nc  =  {id Ax)3.  The  isotropic  turbulence  in  the  periodic  box  has  the 

following  characteristics:  integral  length  scale  L=0.594,  Taylor  micro  length  scale  X 

=0.277,  Kolmogorov  length  scale  r)=0.023,  Kolmogorov  time  scale  tk  =(v/s)1/2  =0.073, 

and  Re  ^  —40. 1 .  The  time  step  used  in  this  study  is  dt=0.01  and  the  ratio  of  total  time 

(after  the  particles  have  reached  dynamic  equilibrium)  to  Kolmogorov  time  is  82  for  most 
of  the  cases  in  the  following  discussions.  The  instantaneous  non-uniform  concentration  is 
calculated  and  recorded  at  every  time  step  after  particles  have  reached  dynamically 
equilibrium.  The  average  concentration  non-uniformity  is  obtained  based  on  the  recorded 
concentration  non-uniformity  by  all  time  steps  and  all  turbulence  realizations. 

Figure  5.7  shows  the  effects  of  inertia  on  the  concentration  non-uniformity  for 
sizeless  particles.  It  can  be  seen  that  the  maximum  value  of  Ap  occurs  around  Pxk  =1.0, 

which  agrees  with  the  results  of  Wang  &  Maxey  (1993).  It  is  noted  that  all  curves  for 
different  inertia  increase  monotonically  as  the  cell  size  decreases,  and  the  slopes  are  all 
close  to  -1  in  the  small  cell  size  limit  except  that  the  curve  for  Pxk  =0.51 1  has  a  slightly 
smaller  slope  at  the  same  region.  Evidently,  the  results  in  figure  5.7  show  the  divergence 
of  Ap  as  cell  size  becomes  small,  so  does  figure  5.6.  The  singular  behavior  of  Ap  is 
unphysical  for  finite-size  particles  because  the  particle  concentration  becomes  either  0  or 
1  per  cell  when  the  cell  size  falls  below  particle  diameter.  Analytically  it  can  be  easily 
shown  from  equations  (5.21)  &  (5.22)  that  Ap  is  equal  to  zero  if  the  0-1  Poisson 
distribution  is  reached. 
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Figure  5.7  The  concentration  non-uniformity  of  sizeless  particles  as  function  of  cell 
size  for  various  inertia  in  isotropic  turbulence. 


Figures  5.8  (a-b)  show  the  non-uniformity  parameter  Ap  for  finite-size  particles 
over  the  whole  range  of  cell  sizes.  The  collection  of  the  maximum  values  of 
concentration  non-uniformity  parameter  Ap  for  finite-size  particles  as  function  of  particle 

inertia  in  isotropic  turbulence  is  shown  in  figure  5.9.  For  Pxk>  0.5,  the  non-uniformity 
parameter  Ap  increases  as  the  inertia  increases.  The  maximum  Ap  is  reached  at 
pxk=0.5.  For  very  large  inertia  Pxk<0.5,  however,  the  non-uniformity  parameter  Ap 
decreases  as  the  inertia  further  increases.  It  is  observed  that  Ap  varies  with  inertia  as 

(Pxk)  for  large  values  of  Pxk.  It  is  also  noted  from  figure  5.9  that  the  non-uniformity 
Ap  at  D/r|=0.957  is  larger  than  that  at  larger  particle  size  D/r|=2.522.  Detailed 
discussions  on  the  relationship  between  the  non-uniformity  parameter  Ap  and  the 
collision  rate  will  be  given  in  Chapter  6. 
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Figure  5.10  compares  Ap  of  small  finite-size  particles  with  that  of  zero  size  for 
Pxk  ~1.  It  can  be  seen  that  there  exists  a  maximum  Ap  for  all  finite-size  particles  for  all 
different  inertia.  Clearly,  the  use  of  sizeless  particles  not  only  results  in  a  singular 
concentration  non-uniformity  in  the  limit  of  small  cell  size  but  also  overestimates  the 
value  Ap  even  at  finite  cell  size.  Figures  5.11  (a)  &  (b)  show  that  the  value  of  the 
maximum  Ap  increases  as  the  particle  size  becomes  small.  A  maximum  value  of  Ap 
always  exists  for  all  finite-size  particles.  It  is  reasonable  to  expect  that  the  maximum  Ap 
will  ultimately  become  infinite  if  one  uses  infinite  small  size  of  particles  such  as  sizeless 
particles.  It  is  seen  that  the  maximum  Ap  is  reached  at  normalized  cell  size  Ax/D=2. 0-2.5 

for  all  particle  sizes  (D/r|=0. 957-2. 522)  except  for  the  smallest  particles  in  figures  5.11 
(a)  &  (b),  where  the  maximum  Ap  is  reached  at  Ax/D=3.5  for  D/r|=0.435. 

The  effects  of  particle  size  on  the  particle  concentration  non-uniformity  may  be 
explained  as  follows.  A  reason  for  particle  size  to  play  an  important  role  in  quantification 
of  the  concentration  non-uniformity  may  be  explained  by  the  mechanism  of  binary 
collisions  between  particles.  Since  particles  are  pairing  to  form  the  non-uniformity  before 
they  collide  each  other  as  turbulent  flow  evolves  in  time,  the  most  meaningful  scaling 
quantity  should  be  able  to  measure  the  degree  of  the  number  of  particles  closely  paired  at 
every  instant.  The  dimensionless  cell  size  Ax/D=2.0~3.5,  with  being  able  to  contain  more 
than  a  pair  of  particles,  is  a  proper  frame  to  capture  this  significant  phenomenon. 
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Figure  5.8    The  concentration  non-uniformity  parameter  Ap  of  finite-size  particles  as 

function  of  cell  size  for  various  inertia  in  isotropic  turbulence,  (a)  D/r)=0.957. 
(b)  D/ti=2.52. 
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Figure  5.9  The  maximum  concentration  non-uniformity  parameter  Ap  of  finite-size 
particles  as  function  of  particle  inertia  in  isotropic  turbulence. 
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Figure  5.10  The  concentration  non-uniformity  parameter  Ap  as  function  of  Ax/r| 
for  small-size  and  zero-size  particles  in  isotropic  turbulence,  px  k = 1 .023 . 
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Figure  5.11  The  concentration  non-uniformity  of  particles  at  a  given  inertia  as  function 
of  Ax/D  cell  size  for  various  sizes  in  isotropic  turbulence,  (a)  Log  scale  in  x-axis. 
(b)  in  normal  scale  in  x-axis. 
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5.4    Summary  and  conclusions 


A  new  method  is  developed  to  calculate  reliably  the  normalized  concentration  non- 

2  ~2 

uniformity  parameter  <n'  >/n  ,  in  which  n  is  the  instantaneous  particle  number 

concentration  and  the  prime  is  the  fluctuation  in  a  small  cube  of  volume  (Ax)3 ,  from  the 

Lagrangian  simulations  of  the  particle  motions  in  turbulence,  and  n  is  the  mean 
concentration  of  particles.  The  inherent  statistical  noise  is  eliminated  from  the  parameter 

2  _2 

<n'  >/n  .  The  effects  of  particle  size  and  inertia  on  the  concentration  non-uniformity 
can  be  quantitatively  examined  using  present  method.  It  is  found  that  the  concentration 
non-uniformity  parameter  will  be  greatly  overestimated  if  the  effect  of  particle  size  is 
neglected.  For  particles  with  zero  size,  the  non-uniformity  varies  with  the  inverse  of  Ax, 
i.e.  Ap  ~  1/Ax  for  small  values  of  Ax  when  particle  response  time  is  comparable  with 
turbulence  Kolmogorov  time  scale.  However,  for  finite-size  particles,  the  concentration 
non-uniformity  parameter  reaches  a  maximum  value  instead  of  unbounded  value  for 
sizeless  particles.  For  finite-size  particles,  the  maximum  non-uniformity  is  reached  at 
Ax/D  =  2.0-3.5  while  it  becomes  zero  when  Ax/D  approaches  to  one.  Particle  size  is  an 
important  length  scale  in  quantifying  particle  concentration  non-uniformity.  The 
maximum  value  of  the  concentration  non-uniformity  is  reached  when  the  ratio  of  particle 
response  time  to  the  Kolmogorov  time  scale  is  close  to  one. 


CHAPTER  6 

COLLISION  RATE  AT  FINITE  PARTICLE  INERTIA  AND  SIZE 

IN  GAUSSIAN  TURBULENCE 

The  effects  of  particle  inertia,  size,  concentration  and  gravity  on  the  collision  kernel  are 
investigated  in  an  isotropic,  Gaussian  turbulence.  An  asymptotic  analysis  is  developed  to 
predict  the  effect  of  the  increased  collision  velocity  on  the  collision  rate  due  to  small 
particle  inertia.  A  unified  expression  of  the  collision  kernel  including  effects  of  the 
increases  in  both  collision  velocity  and  concentration  non-uniformity  is  proposed  for 
small  inertia  particles.  Good  agreement  is  obtained  between  the  simulation  results  and 
the  prediction  by  the  unified  expression  for  particles  with  small  inertia.  The  theoretical 
prediction  of  the  collision  kernel  for  zero  inertia  particles  is  critically  examined  by 
carrying  out  the  direct  numerical  simulations  at  both  a  very  dilute  concentration  with 
small  number  of  particles  and  a  dilute  concentration  with  large  number  of  particles.  It  is 

found  that  the  normalized  collision  kernel  a*n  decreases  as  particle  size  or  mean 
concentration  increases.  In  contrast  to  Saffman  &  Turner's  prediction  for  the  effect  of 
particle  gravity  on  a*! ,  the  gravitational  settling  is  found  to  reduce  a*n  for  monodisperse 
particles.  It  is  found  that  "HARD  SPHERE"  model  suffers  from  multiple  collision 
problem  and  the  collision  kernel  predicted  by  the  "HARD  SPHERE"  model  is  dependent 
on  time  step  size  for  small  inertia  particles.  The  classical  result  of  Smoluchowski  for  the 
collision  rate  of  monodisperse  particles  in  a  laminar  shear  flow  is  shown  to  be  inaccurate 
due  to  the  inclusion  of  the  self-collisions.  The  present  work  extends  Smoluchowski' s 
result  by  excluding  the  self-collisions  in  the  counting  of  collision  pairs. 
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6.1  Introduction 

As  pointed  out  in  Chapter  1,  the  particle  inertia  affects  the  particle  collision  rate 
through  various  mechanisms:  increase  in  the  collision  velocity,  particle  concentration 
non-uniformity,  and  compressibility  of  the  particulate  phase  in  the  continuum  sense. 
Zhou  et  al.  (1997)  made  an  attempt  to  include  all  these  mechanisms  in  a  mechanistic 
analysis  for  particle  collision  rate.  In  particular,  the  effects  of  the  concentration  non- 
uniformity  were  modeled  with  an  adjustable  parameter  that  is  based  on  the  simulation 
results  of  the  collision  using  the  "OVERLAPPING"  scheme.  It  must  be  pointed  out  that 
the  "OVERLAPPING"  scheme  for  particle  detection  scheme  is  equivalent  to  having  a 
zero  particle  size  as  far  as  the  particle  concentration  non-uniformity  is  concerned.  It  has 
clearly  been  demonstrated  in  Chapter  5  that  neglecting  particle  size  leads  to  a  significant 
over-estimation  of  the  concentration  non-uniformity.  Hence,  the  model  of  Zhou  et  al. 
(1997)  contains  a  serious  over-estimation  on  the  effect  of  particle  inertia  through  an  over- 
estimation  of  the  concentration  non-uniformity. 

In  addition  to  particle  inertia  measured  by  Pxk,  the  particle  size  also  plays 
complicated  roles  on  the  collision  rate.  Consider  the  zero  inertia  limit  for  example.  For 
particles  with  small  size,  D<r\,  the  collision  velocity  w  is  linearly  proportional  to  the 

particle  size  D  and  (s  /  v)1/2.    This  is  so  because  for  |Xj  -Xj|~D<r|,  Taylor  series 

expansion  can  be  applied  to  obtain  fluid  velocity  difference,  Au,  between  two  close  points 
separated  by  a  distance  |Xj  -  Xj  |=R  in  space  in  the  form  of 

<|w|>  ~  <|Au  |>  ~  |  ( Xj  -  xj )-Vu  +  ...  |  ~  R(s  /  v)1/2+....  (6.1) 

On  average  the  leading  order  term  for  w  is  thus  proportional  to  (s  /  v)l/2  and  D.  For 
|Xj  -  Xj|  ~  D»r|,  the  velocity  difference  Au  between  two  points  in  space  is  given  by  the 
inertial  subrange  scaling  law  as 

|Au|~  (eR)1/3+...  (6.2) 
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which  is  independent  of  viscosity.  Note  that  the  dependence  of  |Au|  on  R  is  only  to  the 
one-third  power  which  is  quite  less  than  the  linear  behavior  for  R<r|.  Hence,  in  addition 
to  reducing  the  particle  concentration  non-uniformity,  an  increase  in  particle  size  also 
reduce  the  rate  of  increase  in  the  collision  velocity.  When  this  is  coupled  with  particle 
inertia  and  finite  particle  volume  concentration,  the  situation  becomes  much  more 
complicated.  A  detailed  numerical  study  on  the  effects  of  particle  inertia,  particle  size, 
and  particle  concentration  on  the  collision  rate  in  turbulence  is  warranted. 

To  that  end,  an  asymptotic  analysis  is  first  developed  to  predict  the  effect  of  the 
small  particle  inertia  on  the  collision  kernel  due  to  the  increase  in  collision  velocity.  The 
effect  of  particle  size  on  the  collision  kernel  is  critically  examined  by  carrying  out  the 
direct  numerical  simulation  at  a  very  dilute  concentration  for  zero  inertia  particles.  It  is 
found  that  the  collision  kernel  decreases  as  particle  size  increases.  The  collision  kernel  is 
found  to  be  affected  by  particle  mean  concentration  of  the  system.  A  unified  expression 
of  the  collision  kernel  including  effects  of  the  increase  in  collision  velocity  and 
concentration  non-uniformity  is  proposed  for  small  inertia  particles.  The  unified 
expression  is  evaluated  through  direct  numerical  simulations.  The  effect  of  gravitational 
settling  on  the  collision  kernel  is  observed  through  DNS  result  and  explained  by  the 
affected  collision  velocity  and  concentration  non-uniformity.  The  results  of  three 
collision  models:  "THROW  AWAY",  "KEEP  ALL"  and  "HARD  SPHERE"  are 
presented  and  compared.  It  is  found  that  good  agreement  between  the  "THROW 
AWAY"  and  "KEEP  ALL"  models  are  reached  and  verified  by  Smoluchowski's  result  in 
laminar  shear  flow,  while  the  "HARD  SPHERE"  model  gives  quite  different  result  due  to 
the  multiple  collision  mechanism.  It  is  observed  through  DNS  that,  even  at  a  large  inertia 
such  as  pxk=l  and  with  long  time  PT=52,  the  collision  kernel  using  "HARD  SPHERE" 

model  appears  continuously  increasing  with  time,  but  the  results  of  both  the  "THROW 
AWAY"  and  "KEEP  ALL"  models  are  quite  steady  from  the  beginning  and  also  in  a 
good  agreement.  The  Smoluchowski's  collision  rate  for  monodisperse  particles  (1.5)  in  a 
laminar  shear  flow  is  found  to  be  inaccurate  due  to  the  inclusion  of  the  self-collisions.  A 
correction  to  Smoluchowski's  formulation  has  been  made  by  excluding  the  self-collisions 
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in  the  counting  of  collision  pairs.  A  numerical  simulation  for  particle  collisions  in  a 
laminar  shear  flow  at  very  low  concentration  is  carried  out  to  validate  the  present  result. 

6.2  Correction  of  Smoluchowski's  theory 

As  mentioned  in  Chapter  2,  the  collision  rate  and  kernel  may  vary  with  the 
collision  models.  In  order  to  minimize  the  effects  of  the  collision  models  on  the  collision 
rate  and  kernel,  numerical  simulation  of  particle  collisions  has  to  be  carried  out  at  a  very 
dilute  particle  concentration  in  the  system.  Although  equations  (1.3)  and  (1.5)  are  used 
usually  for  systems  involving  a  large  number  of  particles  so  that  the  number 
concentration  in  the  continuum  sense  is  well  defined,  the  calculation  of  the  collision  rate 
given  by  equation  (1.3)  is  not  necessarily  limited  to  the  systems  with  large  number  of 
particles.  There  has  been  recent  development  in  the  numerical  studies  of  particle 
collision  in  fluids  in  which  a  finite  number  of  particle  trajectories,  instead  of  continuous 
particle  cloud,  are  studied.  It  has  come  to  our  attention  that  equation  (1.5)  involves  the 

counting  of  self  collisions  among  monodisperse  particles,  despite  the  factor  of  ^  which 

has  already  been  introduced.  This  issue  can  be  easily  seen  in  the  Lagrangian  framework. 

Suppose  there  are  Nj  particles  of  /th-size  group  and  Nj  particles  of  y'th-size  group 
in  the  fluid  of  volume  V  at  a  given  instant.  The  number  of  collisions  among  /th-size  and 
/th-size  particles  is  related  to  the  number  of  possible  combinations  among  these  two 
groups.  The  number  of  possible  combinations  among  /th-size  and  yth-size  groups  is  NjNj 
so  that  the  collision  rate  can  be  related  to  the  number  of  possible  combinations  as 

nClJ  =  otij  NiNj  /V2.  (6.3) 
However,  considering  the  collisions  among  the  same  size  group,  the  number  of  possible 
combinations  of  particles  is  ^Nj(Nj-l),  not  ^  N? .  Take  Nj=2  for  example.  There  can 
be  only  one  possible  collision  in  the  system.  However,  equation  (1.5)  implies  that  there 
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are  two  possible  collisions.  Based  on  the  number  of  possible  combinations,  the  collision 
rate  among  the  same  size  group  in  the  given  volume  should  be 

nca  -  jNiCNi-UylW3*4  (6-4) 

or 

nca  =  I  ni(ni  -1/V)  i  r(2rj)3  for  n{  >1/V  (6.5) 

where  nj  =  Nj/V  and  nj  =  Nj/V  are  simply  the  average  number  of  particles  per  unit 
volume.  For  a  large  value  of  Nj,  or  njV»l,  equation  (6.5)  recovers  equation  (1.5). 

Using  the  number  of  possible  combinations,  equation  (6.4)  can  also  be 
consistently  derived  from  equation  (6.3)  as  a  special  case.  Considering  Nj=2M  particles 
of  radius  rt  in  a  given  volume  V  and  separating  these  particles  into  two  groups  each 
having  2M~'  particles,  the  number  of  collision  candidates  among  these  two  groups  is 

22(M- 

\  However,  there  are  collisions  among  the  particles  within  each  group.  To  find  the 
number  of  collisions  within  the  group,  all  the  particles  in  the  group  can  be  further  divided 
into  two  sub-groups,  each  with  2M"2  particles  on  such  second  level  sub-group.  The  total 
number  of  the  collision  candidates  in  the  second  level  sub-groups  is 

2x22(M-2)_  This 

process  can  be  continued  until  there  is  only  one  particle  in  each  sub-group  on  the  last 
level  so  that  all  possible  combinations  are  counted.  Summing  up  the  collisions  among  all 
sub-groups  of  different  levels  yields 

nc,,=  ^r(2ri)3  [22(M-1)  +  2x22(M"2)  +  22x22(M'3>+  +  ...  +  2M'2  22  +  2^2°] 

=  lr(2ri)32M(2M-l)/2=  i-Ni(Ni-l)  |r(2rj)3.  (6.6) 

Thus  equation  (6.4)  is  now  consistent  with  equation  (6.3)  and  is  indeed  a  special  case  of 
equation  (6.3). 

To  confirm  the  validity  of  equation  (6.4)  in  the  context  of  particle  collision  in 
fluid  flows,  a  direct  numerical  simulation  is  carried  out  to  examine  the  collisions  among 
the  like  particles.  Fifty  particles  of  radius  rj=0.02  are  randomly  introduced  at  t=0  into  a 
cubical  box  with  length  Lx=2,  Ly=1.8,  and  Lz=2.0  so  that  V=7.2,  with  the  restriction 
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0.1<yp<1.9.  A  uniform  shear  flow  with  a  shear  rate  r=1.0  is  imposed  and  the  flow 
velocity  field  is  given  by  v=(u,  v,  w)  =  (ry,  0,  0).  All  particles  move  with  the  fluid  along 
the  x-direction.  The  motion  of  particles  is  updated  with  a  time  step  At=0.03  and  a  total  of 
three  hundred  time  steps  is  used  in  each  realization.  Periodic  boundary  conditions  are 
employed  along  x-  and  z-directions  so  that  the  particles  that  have  left  the  box  at  the  end 
of  the  nth  time  step  are  reintroduced,  at  (n+1  )th  time  step,  by  image  particles  entering  the 
box  from  the  opposite  side.  In  a  uniform  shear  flow,  no  periodic  condition  can  be 
imposed  along  the  y-direction  and  particles  near  y=0.1  and  y=1.9  have  less  collisions 
because  there  are  no  particles  in  regions  of  y<0.1  and  y>l .9  in  direct  numerical 
simulation.  A  correction  on  the  collision  rate  due  to  this  boundary  effect  (Wang  et  al 
1997)  is  thus  employed  for  the  number  of  collisions 

nc„  =  [nCl,]-    /(!-—  — )  (6.7) 
<~»    l  I,,,  Jsimu  v     16  L 

in  the  shear  flow.  There  is  another  possible  correction  (Wang  et  al  1 997)  due  to  the  finite 
volume  fraction  of  the  particles  in  the  volume.  However,  the  volume  fraction  in  the 
present  case  is  0.0233%  so  that  the  maximum  volume  correction  suggested  is  less  than 
0.1%  according  to  the  correction  scheme  (Wang  et  al  1997);  this  correction  is  neglected 
here.  An  efficient  scheme  (Chen  et  al  1995)  is  used  to  search  for  particle  collisions.  The 
"THROW  AWAY"  collision  model  is  employed.  However,  since  the  collision  kernel  is 
so  low  considered  in  this  case  that  the  probability  of  having  two  collisions  in  the  entire 
300  time  steps  within  one  realization  is  very  small  so  that  the  treatment  of  the  collided 
particles  is  not  a  concern  at  all  in  the  present  case. 

Based  on  equation  (6.4)  and  taking  Nj=50,  it  is  seen  that  the  average  number  of 
collisions  within  0<t<T=9  in  each  realization  is 

nCiiVT=|Ni(Ni-l)  i r(2rl)V1T  =  0.13067.  (6.8) 

In  order  to  achieve  a  reliable  statistics,  10,000  collisions  are  desirable  so  that  a  total  of 
NR=  100,000  realizations  are  used  to  generate  an  estimated  collisions  Nc=13,067  based 
on  equation  (6.8). 
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From  the  direct  numerical  simulation  of  Nr=  100,000  realizations,  with  a  1.31% 
boundary  correction  given  by  equation  (6.7),  the  total  number  of  collisions  is  Nc=13,093. 
This  agrees  with  our  estimate  (Nc=l  3,067)  within  0.2%.  If  equation  (6.10)  is  used,  the 
estimated  number  of  collisions  would  be  Nc-13,333  which  is  nearly  2%  larger  than  the 
simulation  result. 

As  seen  from  equation  (6.4),  if  there  is  more  than  one  collision  in  0<t<T=9,  Nj  on 
the  right  hand  side  will  not  remain  a  constant,  nor  does  nCii .    However,  the  ratio, 

nCii  v  /[  — Nj(Nj-l)]  which  gives  the  collision  kernel  cc^  should  remain  constant.  Based 

on  the  counting  of  the  number  of  collisions  in  each  time  step  from  the  simulation,  the 
collision  kernel  ctjj  can  be  calculated  using 

nc  V 

aii=1  &  .  (6.9) 

On  the  other  hand,  if  equation  (6.10)  is  applied,  one  would  obtain,  numerically,  a 
different  collision  kernel  by 

nc  V 

o.-y-51  •  (6.10) 

-N.N.V"1 
2    1  1 

Figure  6. 1  (a)  shows  the  numerically  obtained  value  of  ctjj  for  0<t<9  .  It  is  seen  that  there 
is  a  significant  noise  due  to  the  small  number  of  collisions  within  each  time  step  despite 
the  fact  that  100,000  realizations  have  been  used.  In  order  to  obtain  reliable  data,  a  time 
average  is  then  employed  for  otjj.  Figure  6. 1  (b)  shows  the  cumulative  time  average  of  otjj 

and  a together  with  the  prediction  by  Smolouchowski  [an  ——  r(2rj)  ].  It  is  clearly  seen 

from  the  time  averaged  result  that  the  present  formulation  using  ^Nj(Nj-l)  gives  a 
collision  kernel  that  is  consistent  with  equation  (6.9)  while  the  original  formulation  using 
-  NjNj  is  consistently  smaller  by  about  2%  comparing  with  that  given  by  equation  (6.9). 
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An  extended  version  is  thus  suggested  for  the  formulation  of  the  population 
balance  equations  in  view  of  the  foregoing  analyses.  The  evolution  equations  of  particle 
population  are 

-ann,[n,  -1/V]+  -a12n,n2  -a.jn.n,-... 

^•a„ni[ni  -a.2nin2  -a22n2[n2  -1/V]+  -a23n2n3-... 

a12n,n2  -a13n,n3  -a23n2n3  -a33n3[n3  -1/ V]+-... 

for  n;V>l,  i=l,  2, 3, ...  (6.11) 

where  [  f  ]  denotes  that  this  term  is  zero  if  f<0.  Mathematically,  n;V  can  be  a  non- 
integer  in  the  population  balance  equations.  The  restriction  on  taking  only  the  non- 
negative  value  of  otjj  n^nj  -  1/V]  ensures  that  there  is  no  generation  in  the  ith  group  due  to 
the  collisions  among  the  like  particles. 

It  needs  to  be  pointed  out  that  the  effect  of  the  correction  in  equation  (6.11)  is 
relatively  minor.  For  a  system  with  large  number  of  particles  (Nj»l),  the  present 
correction  carries  no  tangible  impact  on  the  evolution  of  the  particle  size  distribution. 
However,  the  present  analyses  place  the  evaluation  of  particle  collision  rate,  as  given  by 
equation  (6.5),  on  a  more  consistent  footing. 


dn, 
~dT 

dn2 
dn3 
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Figure  6.1  Collision  kernel  a  {i  against  time  in  laminar  shear  flow.  T=  1.0,  D=0.04. 
(a)  Variation  of  ensemble  average,  (b)  Cumulative  average. 
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6.3  An  Analysis  on  the  Effect  of  Particle  Inertia  on  the  Collision  Rate 
6.3.1  Saffman  &  Turner 's  Analysis 

Saffinan  &  Turner  presented  the  first  theory  for  inertialess  particles  which  follow 
fluid  motion  completely.  The  collision  rate  among  particle  groups  1  and  2  is  evaluated  as 

nCl2=n1n2<JWr<0  -wrdS>  (6.12) 

where  wr  is  the  radial  component  of  the  relative  velocity,  -n(n2  wris  the  particle  flux 
moving  inward  to  the  target  particle,  and  the  integration  is  carried  over  the  entire 
spherical  surface  of  radius  R,  as  shown  in  figure  4. 1 .  Because  of  the  continuity  of  fluid 
flow, 

nc12=^nln2  ^entire  sphere  lwrl<iS> 

=  ^nln2  fentire  sphere   <lWrl>  &  •  (6-13) 

The  integration  in  right  side  of  equation  (6.13)  may  be  expressed  by  equation  (4.2)  for 
isotropic  turbulence,  then  the  collision  rate  yields 

nCl2  =  njn2  2nR2  <|wx|>  .  (6.14) 

Recall  that  x  is  the  local  coordinate  parallel  to  the  line  connecting  the  centers  of  the 
colliding  particles  and  <  >  denotes  ensemble  averaging.  It  should  be  noted  that  for  the 
particles  with  zero  inertia,  the  ensemble  averaged  collision  velocity  <|wx|>  is  expressed 
by  (4.3),  so  that 

nCl2  =  n,n2  R3(^  V2=l.294n,n2  R3£)v2  =  a12  n,n2  (6.15) 

where  a12=1.294R3(-)1/2  .  (6.16) 

v 

For  particles  with  inertia  and  gravity,  Saffman  &  Turner  gave  their  second 
prediction  for  identical  particles 
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a,,=1.67R3(-)1/2  (6.17) 
v 

An  inconsistency  occurs  that  the  constant  1 .67  in  equation  (6.17)  differs  with  the  constant 
1.294  in  equation  (6.16)  as  particle  inertia  goes  to  zero. 


6.3.2  Effect  of  Small  Inertia  on  the  Collision  Kernel 

To  understand  the  effect  of  inertia  on  particle  collisional  rate,  we  consider  an 
asymptotic  theory  for  small  inertia.  It  is  assumed  that  the  particle  size  is  much  less  than 
the  Taylor  micro-length  scale,  and  the  gravitational  force  is  negligible.  As  in  Saffman  & 
Turner  (1956),  we  also  assume  that  turbulence  is  isotropic  and  Gaussian. 

Since  the  collision  velocity  |wx|  is  more  relevant  to  the  collision  kernel  as  seen 
from  equation  (6.13),  we  start  with  |wx|  to  examine  the  effect  of  small  particle  inertia. 
Based  on  the  asymptotic  analysis  of  collision  velocity,  the  mean-squared  collision 
velocity  is  expressed  by  equation  (4.19)  for  small  inertia,  and  the  net  increase  of  the 
mean-squared  collision  velocity  due  to  the  inertia  is  given  by  the  equation  (4.20).  Notice 
that  the  collision  velocity  is  the  relative  velocity  of  two  particles  upon  contact  and  hence 
w x (=  v ix  _  v2x )  *s  a'waYs  a  negative  value  by  projecting  the  relative  velocity  ( v  {  -  v 2 ) 
along  the  radial  direction  ( x  i  -  x2 ), 

wx|X2_X]=±R  <  0  for  all  collisions,  (6.18) 

and 

<wx>|X2_xl=±R  <  0.  (6.19) 

Although  u  and  v  are  Gaussian,  wx  does  not  follow  Gaussian  distribution 
because  of  equations  (6.18)  &  (6.19),  and  it  is  demonstrated  in  the  figures  4.4  (a)  &  (b), 
the  PDF  of  collision  velocity  for  various  inertia.  Clearly,  wx  is  not  Gaussian  mainly 
because  the  curves  of  the  PDF  distributions  are  not  symmetrical  to  the  y-axis  (the  vertical 
line  wx=0).  However,  careful  observation  also  tells  that,  for  small  inertia,  the  PDFs  of 
the  collision  velocity  are  almost  symmetrical  to  the  vertical  line  at  wx=<wx>  and 
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appear  close  to  be  Gaussian  distributions  with  respect  to  wx=<wx>.  It  is  also  noticed 
that  the  collision  velocity  PDFs  of  small  inertia  particles  peak  at  the  same  wx .  That  is  to 
say,  if  the  mean  collision  velocity  <  wx  >  at  zero  inertia  is  taken  out  from  the  raw  w  x ,  the 
remained  is  A  wx  =  wx(P)  -  <wx(P  ->  oo)>,  the  net  increase  of  the  collision  velocity  due 
to  the  inertia,  and  it  is  quite  close  to  symmetrical  and  may  be  approximated  by  Gaussian 
distribution.  Thus  we  have 


<|Awx|>=J-<Aw2  >,  (6.20) 


where  <  Awx  >  is  the  net  increase  of  the  mean-squared  collision  velocity  due  to  the 
inertia  as  defined  by  equation  (4.20)  for  identical  particles.   Substituting  <  Awx  >  by 


equation  (4.20)  into  (6.20), 
-r  -  <( 

7t     p       ,A  Dt       Dt  pz      Dt  '  \'D 


1 2  r  2    .  ,,Dui„    Du2x  „       2    ,Dux.2  R2-,i/2 

<i  a  w  x  i> = j-  [-  -  <(  ulx  -  u2x  x     -     )>+^j<  (-^-r  >-rv 


(6.21) 

Denoting  a,,(p)  as  the  collision  kernel  at  arbitrary  inertia  and  an(p->oo)  as 
the  collision  kernel  at  zero  inertia,  the  net  increase  of  the  collision  kernel  for  identical 
particles  due  to  the  inertia  is  defined  as  5ai{, 

8a,,=  a11(P)-a11(P^-oo) 

=  2tcR2  <|Awx|>,  (6.22) 
and  denoting  the  normalized  collision  kernel  as 

,  (6.23) 

R3(V2 
v 

then  the  net  increase  of  normalized  collision  kernel  can  be  expressed  as  8a*n , 

oan=-^-.  (6.24) 

R3(V2 
v 

Combining  (6.21),  (6.22)  &  (6.24),  we  have 
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*     &      2  Dulx    Du2  2       Du   2  1/2 

(6.25) 

and  the  net  increase  of  collision  rate  for  the  system  of  a  large  number  of  particles  due  to 
particle  inertia  may  be  written  as 


•        1  2    8tv    2    2  ^Dulx  Du2x 

onc—n,  ^—  R  [-  -<(ulx  - u2x )(  —  -—)> 


2     ,Dux,  R^ll/2 


For  non-identical  particles  1  and  2,  the  net  increase  of  the  normalized  collision  kernel  is 
„  .     [87ivr  1     .Dux.2  ,  1     1  ,2     2     .  . /  1  Dulx     1  Du2x 

6aili-[^<(^f)  >(p7-p7)  -  V^-^V^'Y^ 

+  ^^<^)2>i1/2'  (6-27> 

P,P2X2D  Dt 

and  the  net  increase  of  collision  rate  due  to  particle  inertia  may  be  expressed  as 


<.  "  871V  „2r   ,0^2    ,  1      1  .  2 

5nc,2=n1n2)|-R^<(i^)2>(---) 

.2<(Ulx_U2x)(±^.±^        ^_  <(^)2>^  ].«.  (6.2g) 

lx     2xA(3,    Dt     p2    Dt        PiP2      Dt  A.2D 


The  present  prediction,  equation  (6.25),  differs  from  the  second  theory  of  Saffman 
&  Turner  in  two  aspects.  First,  oa^-^O  in  the  limit  p->«>,  which  is  consistent  with  the 
result  of  inertialess  particles,  while  Saffman  &  Turner's  second  prediction  has 
a|j  =1.67  which  is  inconsistent  with  the  first  result.  Second,  the  present  analysis 
includes  two  additional  terms  for  inclusion  of  the  effect  of  the  inertia,  as  seen  from 
equation  (6.25).  The  effect  of  inertia  does  not  vanish  in  (6.25)  while  it  does  in  Saffman  & 
Turner's  prediction  when  P,  =  P2,  shown  in  equation  (6.17).  These  two  terms  would 

have  disappeared  if  — -  was  assumed  to  be  perfectly  correlated  or  XD— »<x>  which  was 
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essentially  assumed  by  Saffinan  &  Turner  (1956,  p.26)  by  stating  that  "the  air  velocities 
near  two  close  particles  are  the  same".  These  two  additional  terms  dominate  the  effect  of 
inertia  on  the  collision  kernel  for  particles  of  nearly  equal  inertia. 

6.3.3  A  Proposed  Expression  on  Collision  Kernel 

As  mentioned  previously,  in  addition  to  the  increased  collision  kernel  due  to  the 
increased  collision  velocity,  the  collision  kernel  is  also  affected  by  the  preferential 
concentration  due  to  the  inertia.  In  Chapter  5,  a  net  increase  of  the  concentration  non- 
uniformity  can  be  quantified  by  the  new  quantity  Ap,  the  equation  (5.22).  To  include  the 

effects  of  both  collision  velocity  and  the  concentration  non-uniformity  on  the  collision 
kernel  due  to  the  inertia,  the  total  collision  rate  for  the  system  of  a  large  number  of 
identical  particles  due  to  the  inertia  may  be  expressed  as 


nCn  =n2(l+Ap)[a;i(P^oo)+5a;,(p)]  R3  A 
-  n  a  n  R  (-) 


(6.29) 


where  n  is  the  mean  concentration,  a  ^(P  ->  oo)  denotes  the  normalized  collision  kernel 
at  zero  inertia,  Sotj^P)  is  the  increased  part  of  normalized  collision  kernel  due  to  the 
finite  inertia,  and       represents  the  overall  normalized  collision  kernel 


at  any  inertia  is  then 


Aa;,=  (l+Ap)[a;i(P^oo)  +  5a;i(p)]-  a;,(P^oo) 
=  Ap[a;i(p^oo)+5a;i(p)]+oc<;i(P). 


(6.31) 


Ill 


For  small  inertia,  5a|1(P)=8a*1,  where  8a*,  is  defined  by  equation  (6.25).  Combining 
the  equations  (5.22),  (6.25)  &  (6.31)  yields  the  final  form  of  the  total  net  increase  of 
normalized  collision  kernel 

Aa;1={[l4i.l]p.[<n,2(P2=co)>-l]p.00}{a;l(^oo) 

n       n  n  n 


fc[     *  <(ulx-u2x)(^-^)>+-^<(^)2>],/2 
s       pR        lx  Dt       Dt  p2X.2D  Dt 


*™  [-  A  <(ulx-u2x)(^-^-)>+  -^<(^)2>]1/2. 
8       pR        lx     2x     Dt       Dt  r     p2X.2D  Dt 

(6.32) 


Equation  (6.32)  is  referred  as  "unified  expression"  of  the  collision  kernel  for  simplicity  in 
the  following  discussions.  Obviously,  Aa*j  is  a  complicated  function  of  the  inertia 
instead  of  constant  or  zero-value  predicted  by  previous  theories. 


6.4    Comparison  of  the  Analysis  with  Simulation  Results 

In  this  section,  only  monodisperse  particles  in  Gaussian  turbulence  are  considered, 
and  "THROW  AWAY"  collision  model  is  employed.  The  variation  of  the  collision 
kernel  with  particle  inertia,  size,  number  concentration  and  gravitational  settling  will  be 
discussed  in  succession. 

6.4. 1   Effect  of  Particle  Inertia 

To  isolate  the  effect  of  the  inertia  on  the  collision  kernel  in  our  simulation  and 
compare  with  the  unified  expression  (6.32),    the  inertialess  part  of  collision  kernel 

a  \ ]  (P  — >  oo)  must  be  taken  out  from  the  raw  collision  kernel  a  \ ,  (P) .  The  inertial  part  of 
collision  kernel  is  then  obtained  as  Aa ,  j = a  l  ]  (p)  -  a  \  l  (p  ->  oo) .  Figures  6.2  (a)  &  (b) 
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show  the  variation  of  the  collision  kernel  Aan  with  the  inertia.  The  contributions  from 
the  preferential  concentration,  Ap,  and  from  the  collision  velocity  increase,  5an,  due  to 

inertia  are  also  plotted.  To  facilitate  a  comparison  between  the  prediction  and  simulation, 
the  unified  prediction  based  on  equation  (6.32)  is  also  plotted  in  figure  6.2  (a),  where  the 
particle  size  is  less  than  Kolmogorov  length  scale,  i.e.  D/r|=0.957.   Figure  6.2  (b)  shows 

similar  results  for  larger  particle  size  at  D/r|=2.522.  The  inertial  increase  &a*n  due  to  the 
increased  collision  velocity  is  calculated  based  on  equation  (6.25)  using  the  multiplication 
factor  obtained  at  zero  inertia.  The  inertial  part  due  to  the  preferential  concentration  is 
obtained  by  taking  the  maximum  of  the  Ap  curves  over  the  whole  range  of  the  cell  sizes 

discussed  in  Chapter  5.  It  is  noted  that  the  application  of  the  multiplication  factor 
obtained  at  zero  inertia  may  be  only  meaningful  for  small  inertia.  The  unified  prediction 

is  then  computed  by  equation  (6.32)  based  on  8a *u  and  Ap  determined.  The  inertial  part 
of  normalized  collision  kernel  Aan  increases  as  the  inertia  increases.  The  inertial 

increase  Aan  is  observed  to  reach  a  maximum  at  (3xk=0.146  and  decreases  thereafter 
for  both  D/r|=0.957  and  D/r|=2.522.  However,  the  maximum  Ap  due  to  inertia  is 
reached  at  (kk=0.5  for  both  D/r|=0.957  and  D/n.=2.522,  as  discussed  in  Chapter  5.  The 

reason  for  the  maximum  Aa|]  is  that  increasing  the  inertia  decorrelates  particles  from 
the  turbulent  flow,  and  both  the  collision  velocity  and  concentration  non-uniformity 
increase,  so  the  collision  kernel  is  enhanced.  However,  the  inertial  increase  will  also 
reduce  the  ability  for  particles  to  pick  up  energy  from  the  turbulence.  Further  increase  of 
the  inertia  will  eventually  disable  particles  to  get  energy  from  the  neighboring  fluid  and 
particles  tend  to  be  randomly  distributed  so  that  the  concentration  non-uniformity  will  be 
reduced  till  disappeared.  Comparison  between  figures  6.2  (a)  &  (b)  also  shows  quite 
different  magnitudes  of  the  increased  collision  kernel  by  different  particle  sizes,  that  is, 
more  enhancement  of  the  collision  kernel  is  observed  as  the  size  becomes  smaller.  It  is 

noted  that  the  predicted  Accn  based  on  the  unified  expression  appears  to  agree  with  the 
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simulation  for  small  inertia  particles.  However,  large  deviation  is  seen  for  Pxk  <  2.  This 
is  expected  since  the  asymptotic  analysis  assumes  small  particle  inertia. 


6.4.2  Effect  of  Particle  Size 


In  order  to  understand  the  effects  of  particle  size  and  concentration,  zero  inertia  is 
first  considered.  Figures  6.3  &  6.4  show  the  collision  kernel  as  a  function  of  the  size  at 
different  range  of  particle  concentration.  Figure  6.3  is  for  very  dilute  particle 
concentration.  Only  500  particles  for  the  size  D/r|=0.435,  300  particles  for  D/r)=0.609, 
200  particles  for  D/n.=0.87  and  50  particles  for  D/r|=0.957,  2.522,  3.478  are  used  in  the 
simulation.  The  reason  for  using  such  dilute  concentration  is  to  avoid  the  dependence  of 
the  collision  kernel  on  the  post-collision  treatment,  as  discussed  in  Chapter  3.  Further 
discussions  on  the  effects  of  the  post-collision  treatment  will  be  given  in  §6.6.  Since  the 
particles  have  zero  inertia  in  figures  6.3  &  6.4,  the  effect  of  particle  size  on  the 
normalized  collision  kernel  can  be  clearly  observed.   The  normalized  collision  kernel 

obtained  by  simulation  is  expected  to  recover  the  theoretical  prediction  (a*] =1.271)  as 

particle  size  approaches  to  the  small  limit  (D/r|«l).  However,  increasing  particle  size 
will  reduce  the  normalized  collision  kernel.  It  appears  difficult  to  investigate  the  collision 
kernel  for  very  small  particles  because  of  enormous  computer  capacity  required,  and  it  is 
also  infeasible  to  investigate  the  collision  kernel  for  large  size  D/r|>5  in  a  "very  dilute 
regime",  because  keeping  "very  dilute  concentration"  may  not  have  enough  particles  to 
generate  reliable  statistics  while  a  given  enough  number  of  particles  may  violate  the 

condition  "very  dilute  concentration".  It  is  seen  from  figure  6.3  that  a\{  approaches 
constant  for  large  particles.  Figure  6.4  shows  the  effect  of  particle  size  on  the  collision 
kernel  for  large  number  of  particles,  where  the  total  particles  are  N=36000  -  500  and  the 
corresponding  sizes  are  D/rpO.217  -  13.36.  Compared  with  figure  6.3,  the  difference  is 
that  the  magnitudes  of  the  collision  kernel  are  smaller  and  continue  reducing  as  the  size 

increases  after  D/n>8.    The  dependence  of  the  normalized  collision  kernel  a*n  on  the 
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size,  shown  in  figures  6.3  &  6.4,  reveals  that  the  simple  cubic  law  (a,,  ~  R3)  can  not 
accurately  describe  the  effect  of  particle  size  on  the  collision  kernel  a  x , . 

6.4.3  Effect  of  Particle  Mean  Concentration 

Based  on  the  results  of  the  effect  of  particle  size  on  the  collision  kernel  in  §6.4.2, 
different  magnitudes  of  the  collision  kernel  observed  between  figure  6.3  and  figure  6.4 

gave  us  a  hint  that  the  average  particle  concentration  may  also  affect  a  ^ .  This  has  been 
confirmed  by  figure  6.5,  where  the  total  number  of  particles  N=500,  1,000,  5,000,  12,000 
and  24,000  with  a  diameter  D/r|=0.435  are  used  for  five  sets  of  direct  numerical 
simulations.  It  clearly  shows  that  the  increase  in  the  particle  concentration  will  reduce 

6.4.4  Effect  of  Particle  Gravitational  Settling 

As  discussed  in  §1.3,  Saffman  &  Turner's  prediction  of  collision  kernel  for 
monodisperse  particles  reduced  to  a  constant  a  ^ =1.67  in  the  presence  of  gravity.  To 
examine  their  results,  a  set  of  numerical  simulations  is  carried  out  with  inclusion  of 
particle  gravity.  In  figure  6.6,  the  effects  of  gravitational  settling  on  the  collision  kernel 

a|j  are  plotted  against  the  settling  velocity  for  the  particles  of  diameter  D/r)=2.52  with 
inertia  at  ptk= 1.023.  As  contrast  with  Saffman  &  Turner's  prediction,  increase  in 
particle  settling  velocity  will  reduce  the  collision  kernel.  The  collision  kernel  drops  down 
quickly  as  settling  velocity  increases  from  0  to  8,  and  further  increase  of  settling  velocity 
causes  the  collision  kernel  to  be  reduced  slowly.  The  reason  for  above  effect  may  be 
explained  as  follows.  The  turbulence  is  consisted  of  numerous  eddies  ranging  from  small 
eddies  to  large  eddies.  Large  eddies  are  responsible  for  advecting  both  the  particles  and 
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small  eddies.  The  relative  velocity  between  nearby  particles  is  mainly  generated  by  the 
small  eddies  which  are  critical  to  the  particle  collisions.  As  the  settling  velocity 
increases,  the  time  period  for  the  particles  moving  through  the  small  eddies  decreases, 
thus  the  effects  of  the  small  eddies  on  the  particles  also  decrease.  Due  to  the  decreasing 
effects  of  the  small  eddies  on  the  particles,  the  relative  velocity  between  nearby  particles 
decreases.  This  has  been  confirmed  by  figure  6.7  (a),  which  shows  that  the  collision 
velocity  indeed  decreases  as  the  settling  velocity  increases.  Decreasing  effects  of  the 
small  eddies  will  also  reduce  the  forces  driving  particles  to  form  preferential 
concentration  seen  in  figure  6.7  (b).  It  clearly  shows  that  the  concentration  non- 
uniformity  varies  with  settling  velocity  in  a  similar  way  as  the  collision  kernel  does. 


6.5  Collision  Rate  in  Polydisperse  Suspension 

In  the  following  simulations,  monosize  particles  are  initially  introduced  into  the 
system,  and  the  "KEEP  ALL"  collision  model  is  applied  to  consider  the  particle 
agglomeration  and  size  evolution.  The  same  turbulence  as  above  is  generated.  The 
number  of  particles  will  decrease  as  the  collision  proceeds  and  large  particles  are  formed. 
The  particle  size  evolution  and  the  variation  of  the  collision  kernel  with  particle  inertia 
and  gravitational  settling  will  be  discussed  as  follows. 

Figures  6.8  (a)  &  (b)  show  the  turbulence  ensemble  average  of  the  collision  kernel 
against  time  for  inertialess  polydisperse  particles.  The  collision  kernel  decreases  at  the 
beginning  and  becomes  stable  after  T/xk>7.0.  The  difference  between  the  simulation 
and  prediction  may  be  due  to  the  effects  of  both  the  mean  particle  concentration  and  the 
size,  as  mentioned  in  §6.4.2  and  §6.4.3.  The  collision  kernel  against  turbulence 
realization  is  given  in  figures  6.9  (a)  &  (b).  Very  smooth  collision  kernels  are  obtained 
by  forty  turbulence  realizations  for  polydisperse  particles  at  the  inertia  Pxk=oo,  1.023, 
0.512  &  0.292. 
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Figure  6.2  The  effect  of  particle  inertia  on  the  increased  collision  kernel  Aa*u  for 
monodisperse  particles  for  (a)  D/r|=0.957.  (b)  D/r)=2.52. 
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Figure  6.3  The  effect  of  particle  size  on  the  collision  kernel  au  for  monodisperse 
particles  at  a  very  dilute  system.  (kk  =  oo  . 
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Figure  6.4  The  effect  of  particle  size  on  the  collision  kernel  a  \ ,  for  monodisperse 
particles  at  'not'  very  dilute  system.  pxk  =  oo  . 
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Figure  6.5  Variation  of  the  collision  kernel  a*u  against  particle  number  concentration 
for  monodisperse  particles.  Pxk=oo,  D/r|=0.435. 
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Figure  6.6  The  effect  of  gravitational  settling  on  the  collision  kernel  a*!  for 
monodisperse  particles.  Pxk=1.023,  D/r|=2.52. 
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Figure  6.7  The  effect  of  gravitational  settling  for  monodisperse  particles.  D/n=2.52, 
Px k =1 .023.  (a)  on  the  collision  velocity,  (b)  on  the  concentration  non-uniformity. 
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Figure  6.8  Collision  kernel  a|,  of  polydisperse  particles  against  time.  pxk  =00, 
D/n.^0.957.  (a)  Variation  of  the  collision  kernel,  (b)  Cumulative  average. 
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Figure  6.9  Collision  kernel  a*!  of  polydisperse  particles  with  various  inertia  against 
turbulence  realization.  D/r|=0.957.  (a)  Variation  of  the  collision  kernel,  (b)  Cumulative 
average  of  the  collision  kernel. 


122 


6.5.1  Particle  Size  Spectrum 

Figure  6.10  shows  the  size  distribution  at  three  time  instants.  The  slope  of  the  size 
evolution  keeps  changing  as  the  collision  proceeds  with  time.  It  appears  that  no  dynamic 
equilibrium  of  particle  size  spectrum  exists  in  present  particle-laiden  system,  while  a 
slope  of -2  was  reported  by  Pearson  et  al  (1984)  and  Balachandar  (1988)  for  the  so-called 
dynamic  equilibrium  of  particle  size  distribution.  It  is  noted  that  a  constant  rate  of 
primary  particles  was  introduced  into  the  system  and  large  particles  were  removed  from 
the  systems  in  Monte  Caro  simulations  of  Pearson  et  al.  and  in  Balachandar,  while 
present  simulation  does  not  add  any  extra  particles  or  remove  any  large  particles.  Figures 
6.11  (a)  &  (b)  show  the  number  evolution  of  both  primary  and  secondary  size  particles, 
where  5,000  particles  of  diameter  D/r|=0.957  are  initially  distributed.  It  is  seen  that  the 
instant  number  of  primary  particles  and  the  generation  of  secondary  size  particles  are 
greatly  affected  by  particle  inertia.  The  population  of  new  particles  grows  quickly  as  the 
inertia  increases. 
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Figure  6.10  The  evolution  of  particle  size  distribution  at  various  instants  for 
polydisperse particles.  D/r|=2.17,  Pxk=1.023. 
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Figure  6.1 1  Inverse  of  total  number  of  particles  against  time  for  polydisperse  particles. 
D/ri=0.957.  (a)  primary  size  ( v0 )  particles,  (b)  secondary  size  (2  v0)  particles. 
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6.5.2  Effect  of  Particle  Inertia 

To  investigate  the  effects  of  the  inertia  on  the  collision  kernel  for  polydisperse 

particles,  the  inertial  part  of  collision  kernel  Acc^a^P)  -  aj,(P  ->  qo)  is  calculated 

first.  Figures  6.12  (a)  &  (b)  show  the  variation  of  the  collision  kernel  Aa||  with  the 

inertia.  The  effects  of  the  increase  in  inertia  on  An,  5a\l  and  the  unified  prediction 

(6.32)  are  also  plotted  in  figures  6.12  (a)  &  (b)  for  comparison.  The  particle  sizes  are 
D/r|=0.957  and  2.52  in  figures  6.12  (a)  &  (b)  respectively.  The  number  of  particles  used 
in  simulation  is  5,000  and  2,000  respectively.  Similar  to  the  monodisperse  suspension, 

the  5a*u  due  to  increased  collision  velocity  is  calculated  based  on  the  prediction  (6.25) 
using  the  multiplication  factor  obtained  by  inertialess  polydisperse  particles,  and  the 
inertial  part  due  to  the  preferential  concentration  is  collected  based  on  the  maximum  of 
the  Ap  curves.  It  should  be  pointed  out  that  the  Ap  curves  for  polydisperse  suspension 

are  obtained  by  calculating  the  concentration  non-uniformity  among  the  primary  size 
particles,  so  that  it  represents  the  pure  contribution  from  the  non-uniform  distribution  of 

primary  particles.  The  unified  prediction  is  obtained  by  8a|j  and  Ap  determined.  It  is 

seen  that  the  inertial  part  of  the  collision  kernel  increases  as  the  inertia  increases.  Similar 

to  monodisperse  particles,  the  inertial  part  Aa*(  is  also  observed  to  reach  a  maximum  at 

pxk=0.146  and  decreases  thereafter  for  both  D/r|=0.957  and  D/n=2.522,  and  the 

maximum  Ap  due  to  inertia  is  also  reached  at  p-ck=0.5  for  both  D/r|=0.957  and 

D/n=2.522.  The  reason  for  the  maximum  Aa^  is  similar  to  the  monodisperse  particles. 
Increasing  the  inertia  decorrelates  particles  from  the  turbulent  flow,  and  both  the  collision 
velocity  and  particle  non-uniformity  increase,  so  the  collision  kernel  is  enhanced. 
However,  increasing  inertia  will  also  reduce  the  ability  for  particles  to  pick  up  energy 
from  the  turbulence.  Further  increase  of  the  inertia  will  eventually  disable  particles  to  get 
energy  from  the  neighboring  fluid  and  particles  tend  to  be  randomly  distributed  so  that  the 
concentration  non-uniformity  will  be  disappeared.  It  is  observed  by  comparison  between 
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figures  6.12  (a)  &  (b)  that  larger  collision  kernel  is  obtained  by  smaller  size  D/r|=0.957. 
Since  all  parameters  used  in  figures  6.12  (a)  &  (b)  are  the  same  except  the  size  of  primary 
particles,  it  is  evident  that  particle  size  has  a  significant  effect  on  the  enhancement  of  the 
inertial  part  Act*! .  It  is  also  noted  that  the  unified  prediction  of  collision  kernel  Actn 
seems  to  be  in  agreement  with  the  simulation  at  small  inertia  for  polydisperse  particles. 
Similar  to  the  monodisperse  particles,  large  deviation  is  observed  for  Pxk  <  2.  This  is 
also  expected  since  the  asymptotic  analysis  assumes  small  particle  inertia.  As  seen  from 
previous  discussions,  the  collision  kernel  is  complicated  by  the  collision  velocity  and 
preferential  concentration  and  possibly  also  by  the  compressibility  of  the  particle  velocity 
field.  Based  on  the  simulation  results  of  both  monodisperse  particles  and  polydisperse 
particles,  however,  it  may  be  concluded  that  the  increases  in  both  collision  velocity  and 
preferential  concentration  due  to  the  inertia  are  the  dominant  factors  enhancing  the 
collision  kernels  for  small  inertial  particles. 

6.5.3  Effect  of  Particle  Gravitational  Settling 

The  effect  of  gravitational  settling  on  the  collision  kernel  a*n  can  be  found  in 

figure  6.13.  The  increase  in  particle  settling  velocity  will  reduce  the  collision  kernel  of 
polydisperse  particles.  The  collision  kernel  varies  rapidly  in  the  region  0<  Vs<8  and  then 
slows  down  and  seems  to  approach  a  constant  at  Vs>8.  The  similar  mechanisms  as 

discussed  in  previous  section  may  be  applied  here  to  explain  the  similar  physics  for 
polydisperse  particles.  Because  of  the  decreasing  effects  of  the  small  eddies  on  the 
particles,  the  collision  velocity  is  expected  to  decrease  as  the  settling  velocity  increases, 
as  shown  in  figure  6.14  (a).  Similarly,  the  concentration  non-uniformity  of  polydisperse 
particles  is  also  expected  to  decrease  as  the  gravitational  settling  is  enhanced,  which  has 
been  consistently  shown  by  the  figure  6.14  (b).  It  is  also  observed  that  small  difference 
between  the  figures  6.14  (a)  &  (b)  and  figures  6.7  (a)  &  (b)  is  consistent  with  the  similar 
comparison  of  the  collision  kernel  between  the  figure  6.13  and  figure  6.6. 
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Figure  6.12  The  effect  of  particle  inertia  on  the  increased  collision  kernel  Act  J,  for 
polydisperse  particles  for  (a)  D/r|=0.957.  (b)  D/r|=2.52. 
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Figure  6.13  The  effect  of  gravitational  settling  on  the  collision  kernel  ct^  for 
polydisperse particles.  Pxk=1.023,  D/r|=2.52. 

6.6  Collision  Rate  in  Hard  Sphere  Suspension  and  Assessment  of  Collision  Models 


As  mentioned  in  chapter  2,  the  "HARD  SPHERE"  model  generates  multiple 
collisions  and  may  not  be  used  for  particles  with  small  inertia.  In  addition  to  the  issue  on 
multiple  collisions,  we  have  found  that  Smoluchowski's  result  in  laminar  shear  flow 
cannot  be  recovered  by  the  "HARD  SPHERE"  collision  model.  Figures  6.15  (a)  &  (b) 
show  the  simulation  results  of  hard  sphere  particles  in  laminar  shear  flow  with  various 
shear  rates  and  inertia.  Figures  6.16  (a)  &  (b)  show  the  comparison  of  three  collision 
models  with  Smoluchowski's  prediction.  While  the  "KEEP  ALL"  and  "THROW 
AWAY"  models  can  give  good  agreement  with  the  Smoluchowski's  result,  the  "HARD 
SPHERE"  model,  however,  underpredicts  the  collision  kernel  by  a  large  percentage. 


Figure  6.14  The  effect  of  gravitational  settling  for  polydisperse  particles.  D/r)=2.52, 
(5xk =1 .023.  (a)  on  the  collision  velocity,  (b)  on  the  concentration  non-uniformity. 


129 


Figures  6.17  (a)  &  (b)  show  the  effect  of  particle  inertia  on  the  collision  kernel 
using  "HARD  SPHERE"  model  in  Gaussian  turbulence.  Maximum  collision  kernel  is 
reached  at  Pxk=2.192  does  not  make  sense,  which  is  mainly  due  to  multiple  collision 
mechanism  as  discussed  above.  The  fact  that  the  collision  kernel  increases  with  time 
even  at  Pxk  = 1.023  may  suggest  that  the  "HARD  SPHERE"  model  may  be  appropriate 
for  simulating  collision  kernel  only  if  Pxk«l.  It  is  seen  from  figures  6.17  (a)  &  (b)  that 
the  collision  kernel  appears  to  reach  steady  value  at  px  k  =0.5 1 2  or  0. 1 46.  Figures  6. 1 8  (a) 
&  (b)  show  the  comparison  of  three  collision  models  for  prediction  of  the  collision  kernel 
a*,  at  Pxk= 1.023.  The  collision  kernels  against  time  can  be  observed  for  the  "KEEP 
ALL"  and  "THROW  AWAY"  models  and  good  agreement  is  also  obtained  between  these 
two  models,  while  the  "HARD  SPHERE"  model  results  in  twice  the  collision  kernel  and 
seems  still  diverging  as  PT=51.  The  comparison  among  three  collision  models  at  a 
larger  inertia,  Pxk  =0.292,  is  given  in  figures  6.19  (a)  &  (b).  The  good  agreement 
between  the  "KEEP  ALL"  and  "THROW  AWAY"  model  is  again  observed.  Although 
the  collision  kernel  obtained  by  using  "HARD  SPHERE"  model  is  steady  at  Pxk  =0.292, 

it  differs  considerably  from  the  results  of  other  two  models. 


6.7  Summary  and  Conclusions 


The  classical  result  of  Smoluchowski  for  the  collision  rate  of  monodisperse 
particles  in  a  laminar  shear  flow  is  shown  to  be  inaccurate  due  to  the  inclusion  of  the  self- 
collisions.  The  present  work  extends  Smoluchowski's  result  by  excluding  the  self- 
collisions  in  the  counting  of  collision  pairs.  A  numerical  simulation  for  particle  collisions 
in  a  laminar  shear  flow  at  very  low  concentration  is  carried  out  to  validate  the  extended 
result  of  Smoluchowski. 
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Figure  6.15  Collision  kernel  aj,  of  hard  sphere  particles  with  various  inertia  against 
time  in  laminar  shear  flow.  D=0.04.  (a)  Variation  of  the  collision  kernel,  (b)  Cumulative 
average. 
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Figure  6.16  The  effects  of  the  collision  schemes  on  the  collision  kernel  a*!  against  time 

in  laminar  shear  flow.  r=1.0,  D=0.04.  (a)  Variation  of  the  collision  kernel,  (b) 
Cumulative  average  of  the  collision  kernel. 
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Figure  6.17  The  effects  of  the  inertia  on  the  collision  kernel  a*,  against  time  for  hard 

sphere  particles  in  Gaussian  turbulence.  D/r|=0.957.  (a)  Variation  of  the  collision  kernel, 
(b)  Cumulative  average  of  the  collision  kernel. 
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Figure  6. 1 8  The  effects  of  the  collision  schemes  on  the  collision  kernel  a  *  (  against  time 
in  Gaussian  turbulence.  D/r|=2.174,  p-ck=1.023.  (a)  Variation  of  the  collision  kernel, 
(b)  Cumulative  average  of  the  collision  kernel. 
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Figure  6.19  The  effects  of  the  collision  schemes  on  the  collision  kernel  of  hard 
sphere  particles  in  Gaussian  turbulence.  D/r|=0.957,  Pxk  =0.292.  (a)  Variation  of  the 
collision  kernel,  (b)  Cumulative  average  of  the  collision  kernel. 
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The  theoretical  prediction  of  the  collision  kernel  for  zero  inertia  particles  is 
critically  examined  by  carrying  out  the  direct  numerical  simulation  at  a  very  dilute 
concentration  with  small  number  of  particles  in  isotropic,  Gaussian  turbulence.  It  is 
found  that  the  simulation  results  of  finite  size  particles  (D/n>0.435)  are  asymptotically 
approaching  to  the  theoretical  prediction.  The  normalized  collision  kernel  a*!  decreases 

as  particle  size  increases,  while  is  a  constant  based  on  all  theoretical  predictions. 
The  direct  numerical  simulations  at  a  dilute  concentration  with  a  large  number  of 
particles  have  also  been  carried  out  for  examining  the  size  effect  on  a*!  in  the 
turbulence.   The  normalized  collision  kernel  a*u  decreases  as  particle  size  increases. 

The  effect  of  particle  mean  concentration  on  a*n  is  examined.   It  is  found  that  a*n 

decreases  as  the  mean  concentration  increases,  while  previous  theories  showed  no  effect 

* 

of  the  mean  concentration  on  a  1 1 . 

The  effect  of  particle  inertia  on  the  collision  kernel  is  investigated  by  both 
analysis  and  numerical  simulations  in  an  isotropic,  Gaussian  turbulence.  A  unified 
expression  of  the  collision  kernel  including  effects  of  the  increased  collision  velocity  and 
concentration  non-uniformity  is  proposed  for  small  inertia.  Good  agreement  is  obtained 
between  the  simulation  results  and  the  unified  expression  for  particles  with  small  inertia. 

In  contrast  to  Saffman  &  Turner's  prediction  for  the  effect  of  particle  gravity  on 
the  collision  kernel.  The  gravitational  settling  is  found  to  reduce  the  collision  kernel 
among  monodisperse  particles.  Sharp  decrease  in  the  collision  kernel  is  observed  when 
the  dimensionless  settling  velocity  increases  from  0  to  8,  and  further  increase  in  the 
settling  velocity  will  also  result  in  the  decrease  of  the  collision  kernel  but  a  slow  rate  of 
decrease. 

The  effect  of  various  collision  models  on  the  collision  kernel  is  also  examined.  It 
is  found  that  the  theoretical  prediction  of  collision  rate  at  small  size  and  zero  inertia  limit 
can  be  recovered  by  both  "THROW  AWAY"  and  "KEEP  ALL"  models,  while  the 
"HARD  SPHERE"  model  cannot  be  used  for  prediction  of  collision  rate  at  small  inertia 
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even  for  laminar  shear  flow  due  to  the  multiple  collision  mechanism.  Reliable 
calculation  for  the  collision  kernel  using  the  "HARD  SPHERE"  model  requires  extremely 
small  time  step  for  small  inertia  particles.  Steady  values  of  collision  kernel  can  be 
reached  only  for  very  large  inertia  using  "HARD  SPHERE"  model.  In  general,  the 
collision  kernel  based  on  the  "HARD  SPHERE"  model  is  larger  than  the  results  obtained 
by  both  the  "THROW  AWAY"  and  "KEEP  ALL"  models  for  particles  with  large  inertia 
in  Gaussian  turbulence. 


CHAPTER  7 
SUMMARY  AND  CONCLUSIONS 


A  new  theoretical  framework  is  developed  to  evaluate  the  collision  rate  and 
collision  velocity  of  very  small  particles  in  general  turbulent  flows.  The  classical  results 
for  the  collision  kernel  by  Smoluchowski  (1917)  for  laminar  shear  flow  and  by  Saffman 
&  Turner  (1956)  for  isotropic  turbulence  are  recovered.  The  present  approach  differs 
significantly  from  that  of  Saffman  &  Turner  in  that  the  ensemble  average  is  taken  after 
the  collision  rate  for  a  given  flow  realization  is  calculated.  This  avoids  the  assumption  of 
isotropy  and  allows  for  the  evaluation  of  the  collision  rate  in  general  turbulent  flows.  The 
collision  velocity  predicted  in  both  laminar  shear  flow  and  isotropic  turbulence  agrees 
well  with  the  results  from  direct  numerical  simulation  of  particle  collisions.  The  present 
theory  is  used  to  evaluate  the  collision  rate  and  velocity  of  small  particles  in  the  rapidly 
sheared  homogeneous  turbulence.  It  is  found  that  the  effect  of  the  turbulence  structure  on 
the  collision  rate  and  velocity  can  be  mainly  captured  by  using  the  mean  turbulent  shear 
rate  (e/v)1/2.  The  combined  effects  of  the  mean  shear  and  the  sheared  turbulence  on  the 
collisions  are  elucidated.  Simple  interpolations  are  developed  to  evaluate  the  collision 
rate  and  velocity  for  arbitrary  mean  flow  shear  rate  and  turbulence  mean  shear  rate. 

As  particle  inertia  increases,  the  collision  velocity  increases  and  particle 
concentration  becomes  increasingly  non-uniform.  These  factors  significantly  affect  the 
particle  collision  rate  in  turbulence.  The  effect  of  inertia  on  the  collision  rate  is  examined 
by  considering  these  two  distinctive  contributions:  the  effect  of  the  increased  collision 
velocity  and  the  effect  of  the  increasing  concentration  non-uniformity. 

Numerical  and  analytical  studies  on  particle  collisions  in  isotropic  and  Gaussian 
turbulence  are  carried  out.  The  effect  of  particle  inertia  on  the  collision  velocity  is 
examined.  An  asymptotic  analysis  is  developed  to  predict  the  effect  of  small  particle 
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inertia  on  the  particle  collision  velocity.  Asymptotic  prediction  for  the  increased  collision 
velocity  agrees  well  with  the  results  of  numerical  simulation.  The  collision  velocity  is 

found  to  contain  terms  that  are  proportional  to  (0x  k )  and  (Pxk )  .  The  effects  of  both 
particle  inertia  and  size  on  the  collision  angle  are  examined  through  numerical 
simulations.  The  collision  angle  depends  strongly  on  the  particle  inertia  and  particle  size. 
The  average  collision  angle  and  the  standard  deviation  of  the  collision  angle  increase  with 
the  increases  of  both  inertia  and  particle  size. 

In  order  to  examine  the  effect  of  concentration  non-uniformity  on  the  collision 
rate  quantitatively,  a  new  method  has  been  developed  to  quantify  the  particle 
concentration  non-uniformity  using  a  concentration  non-uniformity  parameter.  The 
concentration  non-uniformity  parameter  is  the  ratio  of  the  variance  of  particle 
concentration  to  the  square  of  average  particle  concentration  with  inherent  statistical 
noise  removed.  The  effects  of  both  particle  size  and  inertia  on  the  concentration  non- 
uniformity  are  quantitatively  examined.  The  concentration  non-uniformity  will  be  greatly 
overestimated  if  the  effect  of  finite  particle  size  is  not  included.  For  sizeless  particles 
with  particle  response  time  comparable  to  the  turbulence  Kolmogorov  time  scale,  the 
measure  of  concentration  non-uniformity  parameter  is  inversely  proportional  to  the  cell 
length  Ax,  in  which  the  concentration  non-uniformity  is  measured,  for  small  values  of  Ax. 
However,  for  finite-size  particles,  the  concentration  non-uniformity  parameter  reaches  a 
maximum  value  instead  of  unbounded  for  sizeless  particles.  For  finite-size  particles,  the 
maximum  non-uniformity  is  reached  at  Ax/D  =  2.0-3.5  while  it  becomes  zero  when  Ax/D 
is  close  to  one.  Particle  size  is  an  important  quantity  that  affects  particle  concentration 
non-uniformity.  The  maximum  value  of  the  concentration  non-uniformity  is  reached 
when  particle  response  time  is  close  to  the  Kolmogorov  time  scale. 

The  effects  of  the  inertia  on  the  collision  velocity  and  concentration  non- 
uniformity  are  combined  to  predict  the  effect  of  particle  inertia  on  the  collision  kernel.  A 
unified  expression  of  the  collision  kernel  including  both  of  these  is  proposed  for  small 
inertia.  Good  agreement  is  obtained  between  the  simulation  results  and  the  unified 
expression  for  particles  with  small  inertia. 
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The  effect  of  particle  size  on  the  collision  kernel  in  turbulence  is  examined  for 
particles  with  zero  inertia  at  a  very  dilute  concentration  limit.  The  normalized  collision 
kernel  aj^  an  /[R?j(e/ v)1/2]  decreases  as  particle  size  increases.  The  theoretically 
predicted  value  for  a*,  is  only  asymptotically  observed  in  the  small  size  and  very  dilute 
limit.  At  higher  concentration,  is  also  observed  to  decrease  with  increasing  particle 
size.    As  the  mean  concentration  increases,  a*u  decreases  which  suggests  that  the 

theoretically  predicted  value  for  a*!  is  only  valid  in  a  very  dilute  limit. 

The  classical  result  of  Smoluchowski  for  the  collision  rate  of  monodisperse 
particles  in  a  laminar  shear  flow  is  shown  to  be  inaccurate  due  to  the  inclusion  of  the  self- 
collisions.  The  present  work  extends  Smoluchowski's  result  by  excluding  the  self- 
collisions  in  the  counting  of  collision  pairs.  The  extended  result  for  the  collision  rate  is 
validated  by  a  numerical  simulation  for  particle  collisions  in  a  laminar  shear  flow  at  very 
low  concentration. 

Finally,  the  effect  of  various  collision  models  on  the  collision  kernel  is 
investigated.  Theoretical  prediction  for  collision  rate  at  small  size  and  zero  inertia  limit 
can  be  recovered  by  both  "THROW  AWAY"  and  "KEEP  ALL"  models,  while  the 
"HARD  SPHERE"  model  cannot  be  used  for  predicting  the  collision  rate  at  small  inertia 
even  in  laminar  shear  flow  due  to  an  unphysically  large  number  of  multi-collisions 
between  the  same  pair.  Calculation  of  the  collision  kernel  using  the  "HARD  SPHERE" 
model  requires  extremely  small  time  step  for  small  inertia  particles.  Steady  values  of 
collision  kernel  can  be  reached  only  for  very  large  inertia  using  "HARD  SPHERE" 
model.  In  general,  the  collision  kernel  based  on  the  "HARD  SPHERE"  model  is  larger 
than  the  results  obtained  by  both  the  "THROW  AWAY"  and  "KEEP  ALL"  models  for 
particles  with  large  inertia  in  Gaussian  turbulence. 
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